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I. INTRODUCTION 


A radar target, acting as a scatterer of a specified 
incident electromagnetic wave, can be considered as a single 
input, single output, linear time-invariant (LTI) system for 
a fixed field observation point. The target can thus be 
considered as a transfer function with poles and zeros. Baum 
demonstrated at the Air Force Weapons Laboratory that a 
target's induced current response to an incident electro- 
Magnetic wave has identifiable poles determined by the 
composition and structural geometry of the target [1]. lie 
1974, Moffatt and Mains proposed that the target's scattered 
field pole locations are independent of the incident 
electromagnetic excitation, including aspect and polarization 
fe] Morgan has proven theoretically that, for the case of 
a conducting target, the scattering response contains complex 
natural resonances which are independent of the incident 
electromagnetic excitation [3]. By determining the poles of 
a target's response, aspect independent target identification 
can be accomplished through the use of electromagnetic natural 
resonances. 

Although the concept of radar target identification 
through the use of natural resonances was first proposed in 


1974 by Mains and Moffatt [2], only recently have signal 


processing techniques been applied to locate the poles in a 
radar target's response in the presence of noise. Kumaresan- 
Tufts [4] and Cadzow-Solomon (5) have each developed 
algorithms which have proven successful in the presence of 
noise. This thesis develops computer routines based upon 
these two algorithms and examines their respective performance 


and appropriateness using a variety of scattering data. 


A. THE PROBLEM 

Since the performance of signal processing methods varies 
under different conditions, a system employed to identify 
targets would possibly reach a decision based on the combined 
output of several signal processing methods. For example, the 
Kumaresan-Tufts and Cadzow-Solomon methods could be used to 
extract poles from the response of scale model targets. The 
information so gathered could be used to build a data base for 
comparison with data similarly obtained in actual field use. 
The results of this system would serve as one input to a 
larger system. Other methods would provide input to the 
system, such as the K-pulse method of Kennaugh [6] and the 
annihilation filter used by Dunavin [7], Morgan and Dunavin 
{8} and Chen [9]. As the name suggests, an annihilation 
filter annihilates the target's poles. A system using the 
annihilation filter concept would contain many such filters, 


each previously designed to cancel the poles of a specific 


known target. In actual field use, a radar target's response 
would be input into each of the filters, and the target 
selected would be that matching the filter whose output 
exhibits the lowest signal energy. 

A system used to identify radar targets would require the 
following concept of employment. First, information required 
by each of the sub-systems would be obtained for every target 
class of concern. In actual field use, this information would 
be compared against actual radar target responses. The system 
would then determine the identity of the target based on the 


input from each of its sub-systems. 


B. BACKGROUND 

Consider a perfectly conducting target illuminated by an 
electromagnetic field. The current induced on the surface of 
this target at a given point must satisfy the magnetic field 


integral equation (MFIE), [10] 


F(x, t)=2nvi! (E,t)+f [ K(E,F,t)F(r, -1EFL) as (1) 
pv 


where nis an outward unit vector normal to the surface of 
= ; ~1 , 

the object, J is the surface current density, H, is the 

incident magnetic field, and K is a Green's function dyadic. 


The entire equation is most eaSily understood as the sum of 


driven currents and "feedback" currents corresponding to the 


cross-product term and surface integral term respectively. 
The term driven by the magnetic field, 2AxH, forms the physical 
optics portion of the total current. Physical optics 
describes the cross~product term as the induced current 
Without interaction with the rest of the body. The Green's 
function kernel describes the current at a point on the object 
due to the feedback of currents from every other point on the 
object, as previously illuminated by the incident field. The 
current at each point is then summed over the surface of the 
object. Note that the surface integral term is of principal- 
value type; the integral excludes the point r=r . 

Once the incident magnetic field is no longer present, 
the solutions of (1) are considered the natural modes of the 


object. These natural modes are of the form, J,exp(s,) - The 


natural resonance frequencies S, are of the forn, 
S,=90,+j, (2) 


where 0, is the damping rate in Nepers/sec and ®©, 18S the 
frequency in radians/sec. The natural resonances of (2) are 
functions of the structural geometry of the object and are 
independent of the incident magnetic field. To understand 
how these natural resonances are unique to the geometry and 
composition of the object, consider a set of points on the 
object previously illuminated by the incident field, so that 


— | 


H=0. The current at a given point in the set is due to the 


infinite number of feedback currents from every other point 
in the set. Recall that these feedbacks are described by the 
Green's function kernel in the integral term of (1). Since 
the set of points previously illuminated is physically located 
on the same object, the infinite number of paths that connect 
a point with all other points in the set is the same for all 
points in the set. The infinite number of paths are unique 
to the structural geometry of the object and correspond 
exactly to the infinite number of paths taken by currents 
which feedback to a given point via the Green's function 
kernel. Finally, the composition of the target determines the 
surface current density on the object. Although an infinite 
number of resonances existS in any object, only a limited 
number of these will be measurably excited by an incident 
field of finite bandwidth. These resonances described in (2) 
appear as complex conjugate pairs in the left-half portion of 
the s-plane. 

In the far-field, the back-scattered response of a target 


to an incident plane wave is of the form 
A (-rp,t)=ake 2s fede, t-le-F/cyas’ 


where c is the speed of light and p is the unit vector whose 


direction matches that of the plane wave's propagation. 


Equation (3) is the result of integrating the current at 
each point on the target surface for a fixed point in the far- 
field. Recall that the current at each point on the target 
is defined by (1). Thus, the back-scattered far-field can 


be obtained by substituting (1) into (3): 


= " = 
H (-rp,t)=u(t-r/c) { He(-x8- 84 2 H,(-rp,t) exp(s,t) } (4) 
nz 


The currents in (1) produce the field in (4). In fact, each 
term in (4) corresponds to the term in (1) which produced it. 
Specifically, the first term in (4) describes the physical 
optics scattered field generated by the 2AxH current which, 
of course, is the first term in (1). Similarly, the second 
term in (4) is produced by the source-free currents defined 
by the second term in (1). Like the current described in (1), 
the field in (4) is the sum of two terms, a driven term and 
a term containing feedbacks. 

The results of (4) can also be seen as two forms of the 
Singularity Expansion Method (SEM) developed by Baum [1]. As 
shown by Morgan [10], during the early-time portion of the 
target's response, the scattered field is composed of the 
physical optics scattered field and a "Class 2" form of the 
SEM expansion. The class 2 SEM expansion corresponds to the 
second term of (4), wherein the coefficients ~H" are (lame. 


varying as the wave passes over the target, since the currents 


producing this portion of the field are integrated over a 
time-varying surface area. At the instant the wave passes the 
last point of the target, the physical optics field vanishes 
and the remaining term in (4) ais produced by constant 
coefficients H,. The coefficients H, are constant at this 
instant since the surface area in the integral in (3) is now 
constant. This instant also marks the transition of (4) from 
a "class 2" SEM expansion of time-varying coefficients to a 
"class 1" SEM expansion of constant coefficients. The 
scattered field due to a plane wave is therefore composed of 
a physical optics term and a class 2 SEM expansion in the 
early-time, and a simple class 1 expansion in the late-time. 

Actual measurement of the scattered far-zone field would 
be greatly aided by knowledge of the transition time of the 
field from early time to late time. From [10], this 
transition for a monostatic radar would occur at At=T+2(D+d)/c 
seconds after radar turn-on. Here, T is the pulse duration, 
D is the target's dimension along the direction of wave 
propagation, d is the distance between the target and the 
measurement point and c is the speed of light. 

The discussion presented in this section was extracted 
from work done by Morgan in [10]. The reader is referred to 
this work for a more detailed treatment of the material in 


this section. 


C.. #HISTORS 

The results of the previous section form the basis for the 
hypothesis that the natural resonances found in the scattering 
response of a target to an incident electromagnetic wave are 
unique to that target. Additionally, only a finite set of 
these natural resonances are measurably excited by a wave of 
finite bandwidth. In 1974, Moffatt and Mains proposed that 
the extraction of resonances from a target's response to 
electromagnetic excitation could be used for target 
identification. This work related to earlier work in 1965, 
when Kennaugh and Moffatt first developed the concept of a 
radar target as a linear time invariant system. Poles in the 
Zz-plane are directly related to the natural resonances of a 


target 
: (5) 


where sis given by (2) and At is the sampling interval in 
seconds. Hence, pole extraction involves resonance 
identification. The use of pole extraction algorithms is 


discussed in the next chapter. 


II. POLE EXTRACTION ALGORITHMS 


The use of pole extraction algorithms to identify radar 
targets is discussed in this chapter. A brief discussion of 
two methods precedes the in-depth evaluation of the Kumaresan- 
Tufts and Cadzow-Solomon algorithms. The evaluation of the 
latter two algorithms occurs in two stages. First, each 
algorithm will be evaluated in its ability to extract poles 
from data with known poles. Some of the data processed was 
generated at various signal to noise ratios by a computer 
program written by Morgan [11]. Additional data was produced 
by Morgan's time-domain thin wire integral equation computer 
program [12]. In the second stage, a side by side comparison 
is made of poles extracted by each method uSing transient 
scattering measurements for a thin wire and for various model 
aircraft. Comparisons between the two methods are made as the 


aspect of the aircraft is varied. 


A. PREVIOUS WORK 
1. Direct Minimization 
The most direct way to determine the natural 
resonances in a target's response is to minimize the mean- 
Square error between the modeled signal and the received 


Signal. In [10], Morgan determined that the late-time target 


response to a radar could be represented as a sum of damped 


Sinusoids given by 


A al 0, 
aCe Es cos (w,t+6,) ac 
The frequency, Wy and damping rate, icone are the same 


parameters found in the natural resonance defined in (2). 


Phase, 6 and amplitude, A,, are the remaining parameters. 


(tes 
The representation in (6) is the sum of an infinite number of 
resonances. The sampled response to an incident wave of 


finite bandwidth can be modeled as 


a . N O,ndt 
y(not) = y =2Aje cos (W,ndt+6, ) (7) 
{=1 
where At is the sampling interval in seconds. The four 


parameters of (7) must be adjusted to minimize the sampled 


mean-square error signal 


en (y -Y_) (8) 


between the actual discrete sampled received signal Y, and the 
modeled signal Wo The processing required in this 
minimization problem is both inefficient and highly non- 
linear. Nevertheless, Chong used this method to process 
mathematically-generated data down to 15.0 dB signal-to-noise 


(SNR) ratio [13]. 
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2. Prony's Method 

As in direct minimization, Prony's approach. to 
resonance classification focuses on the late-time portion of 
a radar target's response. However, linear processing and 
root solving are used. The late-time response is modeled as 
the output of an LTI system of order oye Each signal received 
at some discrete sample, n, is considered to be the weighted 
sume of K, previous signals. Thus, the finite term 
approximation of the received late-time signal, ys is defined 
by 


Ya Pn (9) 


The z-transform of (9) is 


=} Kp7! 


K K C10} 
Za Z = =a ©) 


b= 


The roots of this polynomial in z are the poles of the system 
model. Therefore, the key to extracting the poles in the 
system's response lies in solving for the coefficients b, of 
oo) . 

A set of Ky+M received signals in M equations (9) can 
be arranged in matrix form as 


oe 
0 b 
SD Kp (11) 


iia 


In Prony's original method, the data matrix Dy is 
exactly determined, and the coefficient vector, b, is solved 
using linear computations. In the presence of noise, Prony 
overdetermines the data matrix by setting M>K,) and solves for 
the coefficient vector by obtaining the least-squares solution 
to the system of equations. 

The Prony method has two major problems. First, the 
poles obtained by the least squares solution to the 
overdetermined matrix may be strongly perturbed by noise [14], 
since noise does not satisfy the causal model of the system. 
Second, the order of the system is generally not known a 
priori. When the estimated order is greater than the actual 
order, poles due to noise are generated. Prony'’s method 
offers no technique for distinguishing between the Signal 
poles and the extra poles caused by overestimation of the 
system's order. If the estimated system order is less than 
the actual order, actual poles are lost and the remaining 


poles are perturbed from their true positions. 


B. KUMARESAN-TUFTS ALGORITHM 

The Kumaresan and Tufts pole extraction algorithm was 
developed by adapting Prony's method to reduce the problems 
addressed in the preceding section. The Kumaresan-Tufts 
algorithm modifies the least-squares Prony method in three 


ways: 


eZ 


1. Processed signals are arranged in a data matrix 
based on a non-casual model of the system. 


ee §6The model of the system is deliberately 
overestimated. 


3. The system of equations determined by the above 
two criteria is solved by using singular value 
decomposition (SVD). 

Kumaresan demonstrates in [15] that the use of singular 
value decomposition tends to force the extra poles of the 
excess-order system inside the unit circle, while the non- 
causal arrangement of the signals tends to force the signal 
poles outside the unit circle. The excess order of the system 
model reduces the effects of noise on the actual poles. Since 
the noise is stationary and stable, it looks the same in 
forward and backward time. 

= Equations 

Recall that in (9), Prony's technique defines the 
received late-time signal as the weighted sum of Kn previous 
Signals, where K, 18 presumed to be the order of the system. 
Kumaresan models the same late-time signal as the weighted sum 
of K) future signals, where K, is greater than the estimated 


order of the system. This non-casual model is given by 


=S pb’ (aU) 
*y PY pegs 
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A system of M such prediction equations can be written in 


matrix form as 


Y e@06.¢6©@ Y e 
Kp by Y 
; ; - 7 yo (13) 
Yuoe Ya ae b, M-1 
Or, in matrix notation, 
D,-b=y (14) 


As in Prony's method, the coefficients 5). are coefficienteson 
a polynomial in z that models the system's late-time response. 
Two simple manipulations of either data matrix leads to the 
relationship between the coefficients of the Prony model and 
the prediction coefficients of the Kumaresan-Tufts model. 
With b,=-1, a prediction coefficient is related to an 


autoregressive coefficient by 


bi= - By (25) 

From the above relationship, it can be shown that the complex 

pole pairs of the causal model are merely conjugate 

reflections across the unit circle of the pole pairs in the 
non-causal model. 

2. Singular Value Decomposition 
The non-causal arrangement of late-time signals in a 
set of system equations, and subsequent processing through 


Singular value decomposition, combine to separate the signal 


14 


and noise into orthogonal spaces. As discussed in the 
preceding paragraph, poles of the non-causal model are 
reflected outside the unit circle. Kumaresan demonstrates in 
[15] that the extra poles of the excess-order system can be 
forced inside the unit circle through the use of SVD. 
Singular value decomposition factors the ! MXK, data 


matrix Dy. into the product of the matrices: 


=Usv! 
lia (16) 


The columns of U (MXM) are eigenvectors of | By), and the 
columns of V (K)XK)) are eigenvectors of DD, hterinees the 
rank of the data matrix, D,. the diagonal matrix yr (MXK, ) 
contains r singular values which are the square roots of the 
nonzero eigenvalues of both D)D,, and Deo By rearranging 
the three matrices in the product, the pseudoinverse of Dy can 


be obtained as 


D*=Vz'U! 
y Cry) 


where 5° is a (K)XM) matrix whose singular values on the 
diagonal are the reciprocals of those in the <= matrix. 


Finally, the coefficient vector b, of minimum Euclidian norn, 


1s given by 
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The coefficient vector b* so obtained is the minimum length 
least-squares solution to (14). In other words, b’ is the 
best possible solution to (14). In the case of noiseless 
data, the extraneous poles generated by the excess-order model 
will always be inside the unit circle when p* is used. This 
result is generally true for noisy data. 
3. Bias Compensation 

Kumaresan and Tufts [4] observed that the addition of 
noise perturbed the singular values of the y matrix of (16). 
If the perturbation of these singular values is not 
compensated, both the signal poles and extraneous poles are 
biased towards the unit circle. Kumaresan and Tufts used a 
compensation method which reduced the bias in their work, but 
did not derive an analytical justification. In [16], Norton 
derived a more valid bias compensation method based on the 
eigenvalue shifting theorem. 

4. Kumaresan and Tufts Compensation 

If the actual order of the system is K,., then the 
first K, singular values of the = matrix in (16) aremaoms 
zero. The remaining Ky-K, singular values are considered 
noise singular values and are zero in the case of noiseless 
data. The addition of noise perturbs the first kK, signal 
Singular values and increases the noise to some non-zero 
value. Kumaresan and Tufts compensated for this increase in 


the singular values due to the noise by subtracting the 
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average of the noise singular values from the signal singular 
values. The noise singular values were then set to zero. 
5. Compensation Based on Eigenvalue Shifting Theorem 

As described in the previous section, the singular 

values of the matrix Dy are the square roots of the 


eigenvalues of eo and DD, . Assume the noisy data matrix 


can be represented by D=S+Ni, where N is composed of the wide- 


sense stationary white noise process v,, given by 


l Kp 
Vy oe -Vuek, 


The expected value of D,Dy, can be obtained by 


Since S is deterministic, E[SS'’]=SS'. Assuming the noise is 
zero mean, the two cross products are zero. Because we assume 
the noise is wide-sense stationary and white, ELNN! J=o2I, 
where of is the noise variance and I is the identity matrix. 
The expected value of 1D Dy thus becomes 


: te c?t (27) 
=oal+ 
E[ D,Dy | : ok 


sey 


Similarly, the expected value of D,D, : the other source of 
Singular values, is 


Tp 1=S'St+o621 (22) 
E({D)D¥1 7 


The assumption in the results of (21) and (22) is that the 
diagonals of E[{N'N)=E({NN'] equals the noise variance o2 .. 
Equations (21) and (22) show that in the mean, the squares of 
the singular values of Dy are increased by the noise variance. 
The results lead to the method of eigenvalue 
compensation recommended by Norton in [16]. Recall from (16) 
that the eigenvalues of Dy are on the diagonal of the 2 
matrix returned by the singular value decomposition of Dy - 
If K, is the actual order of the system, and kK, is the 
estimated order of the system then the remaining Ky—Kp. 
Singular values of the = matrix can be squared and averaged 
to obtain an estimate of the noise variance, of . These noise 
Singular values can then be set to zero. The first K) 
Singular values of the £ matrix are then squared and reduced 
by subtracting the estimate of the noise variance. The square 
root of the difference becomes the new first K, singles 
values of the compensated -2 matrix. Calculations according 
to (17) and (18) can then be carried out in a normal manner 
to obtain poles in the presence of the noise. Eigenvalue 
compensation requires an estimate of the actual order of the 
system. Methods to obtain this estimate set discussed in 


Chapter III. 
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6. Performance 
The Kumaresan-Tufts algorithm was programmed in 
Fortran and tested on various types of data. The program 
appears in Appendix A. 
a. Synthetically Generated Data 
The starting point for evaluating the performance 
of the Kumaresan-Tufts algorithm was with synthetically 
generated data of the form given by (8) and shown here again 


for convenience 


ee OG ndt 
Yaa Oe cos (w,not+6, ) (8) 
r=" 
Again, A,19,,W,,8,, are the amplitude, damping rate, frequency 
and phase of a set of N damped sinusoids. Noisy data was 


created by adding stationary white noise. 
1. Noise Performance 
The algorithm was evaluated at various SNR‘s, 
ranging from 90.0 dB to 7.0 dB. These SNR's are ratios of 
Signal energy to noise energy rather than the ratio of signal- 
to-noise power. Synthetic data so generated more closely 
resembles the exponential decay of signal power typical in 


actual radar measurements. 


ne 


Figure 1 shows the signal produced by two s- 
plane poles at 90.0 dB. Figures 2 through 6 depict the poles 
extracted from this signal at SNR's ranging from 90.0 dB to 
7.0 aB. Obtained poles are shown at their positions within 
the upper right hand quadrant of the unit circle in the z- 
plane. Not shown are conjugates of each pole which are 
located below the real axis outside the figure boundaries. 

Figures 2 through 6 demonstrate outstanding 
performance on noisy data, even at SNR's of 7.0 dB. The 
scaling needed to show a discernible difference between 
results obtained at 30.0 dB and 7.0 dB would necessarily 
exclude one of the poles from the enlarged figure. The 
average distance of the trial poles obtained in the 7.0 dB 
SNR signal from the true poles is on the order of 107%. This 
magnitude corresponds to that of the average estimate of the 
noise variance obtained in successive trials with this signal. 
The correlation between the distance of trial poles from true 
poles and the noise variance estimate was consistently 
observed with each of the different signal-to-noise ratios 
used. Figure 7 depicts the signal of Figure 1 severely 
corrupted by noise having 7.0 dB SNR. 

As discussed previously, the signal-to-noise 
ratio used in the synthetically generated data is the ratio 
of energy. Figure 8 depicts the results of pole extraction 


from the signal shown in Figure 7, but with a late-time 
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Figure 2. Kumaresan-Tufts Poles, Synthetic Data, 90.0 dB SNR 
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Figure 3. Kumaresan-Tufts Poles, Synthetic Data, 30.0 dB SNR 
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Figure 4. Kumaresan-Tufts Poles, Synthetic Data, 20.0 dB SNR 
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Figure 6. Kumaresan-Tufts Poles, Synthetic Data, 7.0 dB SNR 
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beginning ten nanoseconds later. Since the SNR is calculated 
over twenty nanoseconds for both signals, the signal power at 
some later time will clearly be less than the power ten 
nanoseconds earlier. The results in Figure 8 show complete 
breakdown of the algorithm's ability to extract poles. The 
trial poles shown are the poles closest to the true poles, and 
yet they are located at positions whose reflections are inside 
the unit circle where noise poles are typically located. 

The preceding results show outstanding 
accuracy for full-length noisy data but a complete breakdown 
of the algorithm for the same signal with a later transition 
to late-time. These initial observations are supported by 
Similar findings presented in this thesis. 

b. Thin Wire Integral Equation Generated Data 

For simple objects such as a thin wire, the radar 
response of that object can be computed by establishing 
boundary conditions on the object and numerically solving the 
integral equations that describe the surface current. Recall 
the magnetic field integral equation given by (ie 
Simulations produced by Morgan's time-domain thin wire 
integral equation computer program [12] were used to evaluate 
the pole extraction algorithm. The excitation waveform used 
is the double Gaussian pulse depicted in Figure 9. This pulse 
1s a wide Gaussian pulse with a ten percent width of 0.3 


nanoseconds subtracted from a narrow Gaussian pulse with a ten 
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Figure 9. Double Gaussian Pulse 
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nanoseconds subtracted from a narrow Gaussian pulse with a ten 
percent width of 0.15 nanoseconds. 

Figures 10 through 13 depict back scattering 
response of a 0.1 meter length thin-wire, having a radius of 
0.00118 meter, computed at various incident aspects, ranging 
from thirty degrees to ninety degrees. The laboratory 
arrangement for actual measurements simulated by Morgan's 
program is described in [17]}. Ninety degrees represents a 
broadside aspect, while thirty degrees represents the incident 
plane wave having nearly grazing incidence on the wire. The 
poles extracted at each of the four aspect angles are plotted 
in Figure 14. In this figure, and those that follow which 
depict extracted poles, the signal poles lie in or on the unit 
circle, and the noise poles lie outside. 

The results obtained with this rigorous numerical 
computation demonstrate the aspect independence of the 
extracted poles using the Kumaresan-Tufts method. Note that 
only half of the poles were obtained for broadside 
illumination; two even-numbered poles can easily be seen 
outside the unit circle. This results because of the physical 
symmetry of both the wire and the incident field, thus 
precluding excitation of odd-symmetric modal currents and 
their associated natural resonances. 

Figure 15 exemplifies the computed back-scattering 


response of the 0.1 meter thin wire corrupted artificially 


on 
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with noise at a 20.0 dB SNR. Figure 16 shows the poles 
extracted at each of the four angles of incidence used 
previously in Figure 14. Poles of Figure 14 at 90° are now 
missing in Figure 16, and only the first three low frequency 
poles are tightly grouped. The loss of high frequency poles 
1s expected because these have the highest damping and thus 
lose their energy at the fastest rate. Further comparison 
between results computed at 20.0 dB SNR and infinite SNR are 
offered, angle by angle, in Figures 17 through 20. 

One additional test of the computed thin wire 
scattering was conducted at a 7.0 dB SNR. The corrupted 
waveforms are exemplified by Figure 21; the extracted poles 
are shown in Figure 22. The number of poles obtained has 
decreased with respect to the number obtained at 20.0 dB SNR. 
The grouping of the clusters has also expanded. Angle by 
angle comparisons are again offered in Figures 23 through 26. 

c. Scale Model Measurements 

The transient scattering measurements of scale 
models used for evaluation in this section were made by Walsh 
using the anechoic chamber of the Transient Electromagnetic 
Scattering Laboratory at the Naval Postgraduate School. The 
entire measurement process and laboratory setup are described 


maedetail in ({i7). 
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Figure 20. Integral Equation Thin Wire Comparison, Noiseless 


vs. 20.0 GB SNR, 90 Degree Aspect 
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Figure 23. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 30 Degree Aspect 
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Figure 24. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 GB SNR, 45 Degree Aspect 
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Figure 25. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 GB SNR, 60 Degree Aspect 
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Figure 26. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 AB SNR, 90 Degree Aspect 
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1. Wire Targets 

The thin wire measurements were obtained from 
the scattering response of a 0.1 meter length thin wire having 
radius 0.00118 meter. Recall that these are the same 
dimensions as the wire whose computed response was processed 
in the previous section. The measurements at each of four 
incident aspects are shown in Figures 27 through 30. 

The poles extracted from the four measurements 
are depicted in Figure 31. As before in the computed noisy 
data, tight clusters occur only at the lowest frequencies. 
The poles in these tight clusters are those which are 
measurably present at various aspects. The poles extracted 
at higher frequencies are those which possessed sufficient 
measurable energy at the given aspect. Figure 32 depicts the 
comparison between poles extracted from the measured and 
computed signals. Again, the closest agreement between the 
two sets of poles occurs at the lowest frequences. 

2. Aircraft Models 
Plastic 1/72 scale aircraft models, coated with 
Silver, were used for transient scattering measurements. 
Representative scattering signatures of two aircraft targets, 
measured at six different aspects, are shown in Figures 33 
Enrough 36. 
The results of pole extraction in target 1 are 


shown for a total of six different aspects in Figures 37 and 


50 


x02 30 degree backscattering 


y(t) (volts) 





0 Z 4 


Baguises 27, 


6 8 10 12 14 16 18 20 


Time (nanoseconds) 


Measured Thin Wire Scattering, 30 Degree Aspect 


om 


y(t) (volts) 


Mcnre 28 . 


Zz 4 6 8 10 2 14 16 


45 degree backscattering 





18 20 


Time (nanoseconds) 


Measured Thin Wire Scattering, 45 Degree Aspect 


a2 


x1 03 60 degree backscattering 


y(t) (volts) 





8 10 2 14 16 


18 20 


Time (nanoseconds) 


Figure 29. Measured Thin Wire Scattering, 60 Degree Aspect 


53 


x1 0-3 


90 degree backscattering 


y(t) (volts) 





4 | 
0 2 4 6 8 10 12 14 16 18 20 
Time (nanoseconds) 
Figure 30. Measured Thin Wire Scattering, 90 Degree Aspect 
C 


54 


Imaginary z 


Figure 31. 





Extracted poles 


° ° 


+ 30 degrees 


x 45 degrees 
* 60 degrees 
o 90 degrees 


Real z 


Kumaresan-Tufts Poles, Measured Thin Wire 


32 


Extracted poles 


+ computed > 


oO measured 


Imaginary z 





Figure 32. Thin Wire Comparison, Measured vs. Integral 
Equation 
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Bo. The poles extracted at all six aspects are shown in 


Figure 39. Only one clearly discernible cluster is present 
in each of the three figures. At higher frequencies, no 
useful information is imparted by the data. Results of 


Similar, though slightly improved quality, were obtained from 
target 2. These results are presented in Figures 40 through 
42 in the format of Figures 37 through 39 respectively. 

Although the Kumaresan-Tufts algorithm is 
capable of extracting low frequency poles acceptably, the 
inconsistent results at higher frequences reveals the inherent 
weakness in an algorithm capable of processing only the late- 
time portion of a target's radar response. 

A side-by-side comparison of poles obtained from 
both aircraft by both the Kumaresan-Tufts method and the 
Cadzow-Solomon method is presented at the end of the chapter 


to illustrate the gains afforded by processing the early—time. 


C. CADZOW-SOLOMON ALGORITHM 

Recall from the results depicted in Figure 8 that a late 
transition to late-time, and the consequent reduction of 
Signal power, caused complete breakdown of the Kumerasan-Tufts 
algorithm. The Cadzow-Solomon algorithm addresses this 
shortcoming by processing the signal at the instantaneous 
onset of early-time. Thus, the Cadzow-Solomon algorithm is 


capable of processing the earliest response of a target to 
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Target 2 


electromagnetic excitation where the response has the greatest 
magnitude. 
1. Applicability 
The early-time portion of a target's scattered field 
occurs as long as there 18S a driven portion of the total 
field. Once the field no longer contains a scattered response 
due, in part, to the incident excitation at points on the 
object, early-time ceases and late—time begins. Hence, the 
Cadzow-Solomon models both the system's input and output, and 
equivalently, the poles and zeros of the system transfer 
Pwnetion. 
2. Equations 
The Cadzow-Solomon algorithm extends the  auto- 
regressive equation (9), used in Prony's method, to the more 
general autoregressive moving average (ARMA) equation 
Kp Ky 


y= by +2 a\X,_ 23) 
n T=! n-} 1=0 


where the second summation term models the excitation to the 


system. 
A set of M such equations in matrix form is given by 
x bk, 
Ms, meee Yen noses ky , 
e e e oy Ky 
2a = (24) 
: : ens “a are : 
ai YK eM-2 M-1 Ky*M-1 Ky Yq eM eI 
ag 


68 


As in the Kumaresan-Tufts method, M is selected to be greater 
than the column dimension of the data matrix which is 
Kyt+Ky+1 . 

3. Excess Poles and Noise Removal 

The cCadzow-Solomon method used in this thesis is a 
modification which incorporates the non-causal arrangement of 
the system equations used by Kumaresan-Tufts. This 
modification was first discussed by Norton in [16]. The 
Kumaresan approach of overestimating the system order can be 
used as before in a non-causal model to constrain the noise 
poles inside the unit circle, while SVD forces the signal 
poles outside the unit circle. 

Since the input waveform is known, its order can be 
almost exactly determined. In all the work of this thesis, 
the input waveform used is the double GaussSian depicted in 
Figure 14. Approximately 25 samples defining this pulse of 
0.5 nanoseconds duration makes K, equal 25 in equation (23). 
Since the input is causal, the signal zeros fall inside the 
unit circle where they cannot be easily segregated from 
Similarly located noise poles. However, the signal zeros 
impart no information about the target and need not be 
extracted. The inclusion of the input in the data matrix is 
nevertheless vital to the model of the system and the accurate 


determination of the signal poles. 


69 


The ARMA equation of (23) can be modified to obtain 
Kp Ky 


Ye Oe + 2 8X1 (25) 


Kytn-itl -=0 


The recursive portion of (25) is now in a non-causal form 
Similar to expression (12). A set of M such equations in 


Matrix form is given by 


De 

y Oe eas v4 x e¢ @ > e u 

Kytl KytKp 0 Kw ve 

e e e e b { ~D 

e e - e "ay = : (26) 
Yan’ YK ekp eM Ryn 0 XK ee 2 Yq M-I 

Oe, in matrix notation 
oe el : =rD-: 7) 

[Py JL -B- J=y wnere © [Pye J=[DyDx] (27) 

4. Singular Value Decomposition 


Like the system equations of the Kumaresan-Tufts 
model, the system equations in (26) are processed using 
Singular value decomposition. The coefficient vector is again 
the minimum-norm solution, which constrains the extraneous 
poles and extraneous zeros to be inside the unit circle. 

5. Bias Compensation in the Cadzow-Solomon Formulation 

By compensating the eigenvalues of the =f: matrix in 
(16), the performance of the Kumaresan-Tufts algorithm is 
Significantly improved in the presence of noise. Cadzow- 


Solomon have shown [5] that if the actual orders Ky and Ky 
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are overestimated to be Kp and Kee min (Ky-Kp » Ky-Ky) Singular 
values are zero in noiseless data. Since the input data is 
known, the eigenvalues of the data matrix may be compensated 
in the same manner as in the Kumaresan-Tufts algorithm for 
noiseless data. 

To understand the compensation required in noisy data, 
an analysis of additive noise is required. As given by Norton 
[16], if the input data noise is w, and the output data noise 


is Vio the data matrix may be modeled as 


eee eal xe. 


(28) 
where 
[Nyx J=LNytNy] (29) 
and 
: K Vieee Vy 
: : = 2 ° (30) 
N=]: : i= : ’ 
Wye : - Wark) Vu aaa *VeeKy 
The expected value of Do is then 
T 
E(D,,Dj,]=Sy.Sy,tE (Ny Nyx] (31) 


If the input and output noise variances are not equal, the 
eigenvalue shifting theorem used in Kumaresan-Tufts cannot be 
used to analytically predict the requisite eigenvalue 
compensation of Delo Me: Nevertheless, when the input and 
output variances were assumed equal, and eigenvalue 


compensation similar to that used in Kumaresan-Tufts was 


Tie 


performed, the results were consistently superior to those 
obtained without compensation. Therefore, the results of 
Cadzow-Solomon signal processing presented in this thesis were 
obtained using eigenvalue compensation and the assumption of 
equal noise variance. 
6. Performance 
The Cadzow-Solomon algorithm was programmed in Fortran 
and tested on the same data used for evaluating the Kumaresan- 
Tufts algorithm. Note that the Cadzow-Solomon algorithm 
can use the early-time portion of the data that the Kumaresan- 
Tufts algorithm can not use. The program appears in 
Appendix B. 
a. Synthetically Generated Data 
The starting point for evaluating the performance 
of the Cadzow-Solomon algorithm was with synthetically 
generated data of the form given by (8) plus the addition of 
input data required to model early time data. 
1. Noise Performance 
The algorithm was evaluated at various signal- 
to-noise ratios, ranging from 90.0 dB to 7.0 dB. Figure 43 
shows the signal produced by two s-plane poles at 90.0 dB, 
with a late-time beginning at 10.0 nanoseconds. Figures 44 
through 48 depict the poles extracted from this signal at the 


different signal-to-noise ratios. 
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Figure 44. Cadadzow-Solomon Poles, Synthetic Data, 90.0 GB SNR 
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Figure 45. Cadzow-Solomon Poles, Synthetic Data, 30.0 dB SNR 
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The figures chart the steady degradation of 
the algorithm's performance with the increase of noise. At 
30.0 dB, the location of the low frequency pole is already 
slightly displaced. More significant is the location of one 
of the extracted poles in the noise signal space. At 20.0 GB, 
the low frequency pole is located in some trials on the real 
axis. At 10.0 dB, all the extractions are located on the real 
axis and at 7.0 dB their locations there are dispersed. The 
extraction of the higher frequency pole is 
uncharacteristically more accurate than that of the low 
frequency pole. Even at 7.0 dB, the high frequency pole is 
located with excellent accuracy. The location of the low 
frequency pole near the real axis was chosen deliberately to 
illustrate the difficulty in resolving the slight frequency 
difference between the true pole and a noise pole located on 
the real axis. Also, fewer points were processed using the 
Cadzow-Solomon method than were processed using the Kumaresan- 
Tufts method, since the largest data matrix allowed by the 
programs in Appendices A and B contain fewer data points in 
the Cadzow-Solomon data matrix than in the Kumaresan-Tufts 
data matrix. The results demonstrate the need to process a 
substantial number of points in order to accurately extract 


low frequency poles. 


TS, 


b. Thin Wire Integral Equation Generated Data 

The performance of the Cadzow-Solomon algorithm 
was evaluated using the same set of data tested by the 
Kumaresan-Tufts algorithm. The results are presented in 
Figure 49. Tight clusters appear at frequencies higher than 
those obtained with the Kumaresan-Tufts algorithm. Figure 50 
depicts the poles extracted from the same signal at a 20.0 dB 
SNR. The clustering at this SNR is comparable to the results 
obtained by the Kumaresan-Tufts method with the noiseless 
data. Further angle-by-angle comparisons of the _ poles 
extracted from the noiseless data and the 20.0 dB data are 
depicted in Figures 51 through 54. Note the small number of 
poles in Figure 54 due to the unexcited odd-symmetric poles 
at 90° aspect. 

One further test was conducted on computed data 
at a 7.0 dB SNR. The results are depicted in Figure 55. Even 
at 7.0 dB, discernible clusters are present. Angle-by-angle 
comparisons of the poles obtained in 7.0 dB data and those 
obtained in noiseless data are presented in Figures 56 through 
59. 

c. Scale Models 

The same scale models used to evaluate the 

Kumaresan-Tufts algorithm were used to evaluate the Cadzow- 


Solomon algorithm. 
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Figure 49. Cadzow-Solomon Poles, Noiseless Thin Wire Data 
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Figure 51. Integral Equation Thin Wire Comparison, Noiseless 
vS. 20.0 dB SNR, 30 Degree Aspect 
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Figure 53. Integral Equation Thin Wire Comparison, Noiseless. 
vs. 20.0 QB SNR, 60 Degree Aspect 


Se 


Extracted poles 





N 
e | 
& + noiseless 
= o 20 dB SNR 
Real z 
Figure 54. Integral Equation Thin Wire Comparison, Noiseless 


vs. 20.0 dB SNR, 90 Degree Aspect 


86 


Imaginary z 


Piguise 55. 


Extracted poles 


x 45 degrees 
Pou deprees 


o 90 degrees 


Cadzow-Solomon Poles, 


87 





v. OnaB ONK 


Extracted poles 





N 
E 
2 . noiseless 
Z o07dB SNR 
= Real z 
Figure 56 Integral Equation Thin Wire Comparison, Noiseless 


vs. 7.0 @B SNR, 30 Degree Aspect 


88 


Extracted poles 


we ef HOISeeSS. ow 


0 7dBSNR 


Imaginary z 





Real z 


Figure 57. Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 dB SNR, 45 Degree Aspect 
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Figure 59 Integral Equation Thin Wire Comparison, Noiseless 
vs. 7.0 @B SNR, 90 Degree Aspect 
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1. Wire Targets 

Figure 60 depicts the poles extracted from 
measurements of a 0.1 meter wire. Three tight clusters appear 
at the lowest frequencies and at the highest frequencies. The 
poles in between can not be easily discriminated. The 
dispersion of these poles is apparently due to the aspect 
dependence of their measurable power. In other words, these 
poles are excited more at some aspects then at others. 

Figure 61 depicts the comparison between poles 
extracted from computed data and measured data. As in Figure 
60, close agreement exists at the highest and lowest 
frequencies. The results are much more favorable than those 
Similarly obtained by the Kumaresan-Tufts algorithm. 

2. Model Aircraft 

Figures 62 through 64 depict poles extracted 
from aircraft target 1. As in the Kumaresan-Tufts testing, 
the Cadzow-Solomon testing was conducted at six different 
aspects. Results for target 2 are depicted in Figures 65 
through 67. The results of both targets show clearly defined 
clusters. The first two clusters of target 2 are 
exceptionally tight. However, the mid-frequency clusters of 
target 2 are not as clearly formed as those of target 1. 

Comparisons of poles obtained with each method 
for target 1 and 2 are depicted in Figure 68 and 69 


respectively. These two figures graphically depict the clear 
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Figure 61. Thin Wire Comparison, Measured vs. Integral 
Equation 
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Figure 62. 
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Figure 63. 
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Figure 64. Cadzow-Solomon Poles hareectel, All Six Aspects 
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Figure 66. 
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Figure 68. Pole Comparisons, Target 1, All Six Targets 
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Superiority of the Cadzow-Solomon algorithm over’ the 
Kumaresan-Tufts algorithm. 

[inieerder te Obtain van Initial indication of 
the possibility for target classification through pole 
extraction, nose-on measurements of two additional aircraft 
models were made, processed and compared with the results of 
targets 1 and 2. The nose-on measurements of targets 3 and 
4 appear in Figures 70 and 71 respectively. A comparison plot 
of poles extracted from each of the four targets is depicted 
in Figure 72. Each of the four aircraft measured are fighters 
of similar size and shape (see Table 1). The poles for each 
target are sufficiently different in this single measurement 
to identify each aircraft individually. However, some of the 
poles are arranged in clusters which appear with a harmonic 
pattern similar to that obtained for either of the first two 
aircraft at various aspects. In order to more fully assess 
the target classification capability of pole extraction, 
several measurements should be made of a given aircraft model. 
A plot of the poles extracted from each of these measurements 
would form clusters at the locations of the true poles. The 
centroid of each of these clusters would then be compared 
against the centroid poles similarly obtained from other 
aircraft. Although several poles of different aircraft might 
be similar, the set of poles belonging to an aircraft could 


form the basis for classification if that set was unique among 


10s 


x 


nN 
ina 


1S) 


0.5 


y(t) (volts) 


103 Nose on 


ee ee ee ee ee ee er ee ee ee ee 
‘ 
eee mee eenee eee. ee ee ee ee ee ee ee ee eee eee ee ee ee eee eee ee ey ee ee ee ee 2 


eTrTerererere ry ee OCC CSCO Ome eH HH Oe seeee somes 


PC eeccccctocccee ng he rhP Brest cestecomecesess a8 


Fae Comme sees ee: peeesns seers rese segs eee sesesssestseogeeesees «+28 


Pee ee meet een cee ae FE Be cece eee eee Oe we OOO IE OOH BH eee HMO HH HOES HEE HEH EOE EEE HH HEH SESE SOH HSER OS ESOS SHOES SEH ESE SENOS SS EHOS EHO EH HHH HHH ESEOES 


SOOO OTOH Hoe oOo HH 0 OOo OOOO ESO HM OE OEE OEOE HEH HHO e AMES EHH H HHH SHS OEe MEE SE ESE SESH COOH OT OE HEE SEER EEE HSH PHOS EH SH SH SEH OH ESE H ORE HE SEH SESES SESH SH SE SEE SE SEEESSSHEHH SEH OH EHS EOOR OE 


2 6 8 10 iZ 14 16 18 20 


Time (nanoseconds) 


Figure 70. Target 3 Scattering, Nose-on 


104 


y(t) (volts) 


ee ee ee 


x1 0:3 


oe oes nee = 


Figure 71. 


Nose-on 


6 8 10 12 14 16 


= Time (nanoseconds) 


Target 4 Scattering, Nose-on 


105 





Extracted poles 


x target 2 


Imaginary z 


* target 3 


O target 4 





Figure 72. Cadzow- Solomon Pole Comparisons, 4 Targets, Nose-on 
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the sets belonging to all other measured aircraft. The results 
in Figure 72 demonstrate the possibility of using the Cadzow- 
Solomon pole extraction algorithm to aid in the classification 
of aircraft, perhaps by use of the extracted poles in 


constructing annihilation filters. 


TABLE 1. FULL SIZE DIMENSIONS OF TARGETS RECORDED 


Target number 1 2 3) 4 


Overall length 12.20 als) 210). 16.94 16.00 
(meters) 


Overall height Sy RSS! 5.09 pg ul 4.80 
(meters) 


Wingspan 10.96 10.00 P11. 43 oS 
(meters) 


Tailplane span Unknown 51.55 6.92 5.72 
(meters) 
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III. SUMMARIES AND CONCLUSIONS 


In this chapter, a step-by-step guide through each 
algorithm is presented. At each step, techniques and lessons 
learned are discussed together with general observations. 


Conclusions are presented at the end of the chapter. 


A. KUMARESAN-TUFTS 

The first step in processing a signal with the Kumaresan- 
Tufts algorithm is to determine the beginning of early-time. 
The objective is to pick the earliest possible starting point 
without entering into the latter part of early-time. If the 
starting point for processing is improperly chosen to include 
the early-time, the results will be completely unreliable 
Since the signal no longer satisfies the late time model. If 
the starting point is chosen too late, the signal may not be 
sufficiently strong in the presence of measurement noise. 
Since the signal is the sum of exponentially damped sinusoids, 
the optimum starting point is at the precise instant of 
transaction into late-time. The key to determining the 
beginning of late-time is in determining the beginning of 
early time. Determining the first response of the target to 
excitation cannot usually be done by a simple visual 


inspection of measurement data. Unless the exact distance to 
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the target is known, the most accurate method attempted by the 
author for determining the beginning of early—time is to 
process the signal using the Cadzow-Solomon algorithm. This 
is discussed in the next section. However, the reliance of 
the Kumaresan-Tufts algorithm on information provided by the 
Cadzow-Solomon algorithm is an obvious disadvantage of the 
former method. 

Once the starting point for processing has been selected, 
the next step is to determine the dimensions of the data 
matrix and, consequently, the number of points in the signal 
to be processed. In trials conducted on noiseless synthetic 
data, the accuracy of pole extraction increased steadily with 
the increase in the data matrix dimensions. These trials were 
conducted up to the limit of the array dimensions defined in 
the computer program of Appendix A. The number of points 
processed in measurement data should be as large as possible, 
while still meeting the following two constraints. First, 
incorporate as many cycles of the data as possible. Usually, 
visual inspection of the data reveals a repeating pattern 
which should be entirely incorporated into the window of 
points to be processed. When only portions of these patterns 
are selected, a disproportionate weighting tends to be placed 
on certain poles. Second, signal portions late in the 
response which are no longer distinguishable in the presence 


of measurement noise should not be selected. 


EO) 


The final step involves determining the number of true 
poles in the system. The following approach has proven to be 
the most successful. First, process the signal without any 
eigenvalue compensation to establish an upper bound on the 
order of the system. In most cases, the number of poles 
outside the unit circle will be less than the overestimated 
order of the system. If not, increase the row dimension of 
the data matrix in order to increase the estimated order of 
the system, and repeat. When the number of poles is less than 
the estimated order of the system, then one should gradually 
increase the number of eigenvalues compensated in successive 
trials, while closely observing the effects induced on the 
poles outside the unit circle. As the number of eigenvalues 
compensated is steadily increased, noise poles and weak signal 
poles will move inside the unit circle. The programs in 
Appendix A and B allow the user to compare the results of 
successive trials, by generating overlays for each plot. If 
N poles are in the signal space, at least the first N 
eigenvalues must not be compensated, or true poles may be 
lost. As the actual order of the system is approached by 
compensation, the user will notice an orderly, even 
arrangement assumed by the noise poles. If certain poles 
still remain suspect after compensation, vary slightly the 


other parameters, such as the starting point and the 
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dimensions of the data matrix. Generally, only true signal 
poles will repeatedly assert themselves under’ varying 


parameters. 


B. CADZOW-SOLOMON 

The techniques and general observations offered in the 
preceding section apply equally to the Cadzow-Solomon 
algorithm. An important consideration in this method, not 
discussed above, is the selection of the beginning of early- 
time. Candidates for a starting point are usually at or near 


zero crossings within approximately thirty points of the 


object's £iyvst definite response to electromagnetic 
excitation. Begin processing at the chosen point while 
varying parameters in successive trials. Select the point 


whose successive results are the most consistent under varying 
parameters. 

The selection of the starting point for beginning of 
early-time can be very critical. For example, not a single 
pole could be extracted in one trial wherein the starting 
point occurred only ten points after the actual starting 
point. Additionally, in most cases observed, the late-time 
start given by the selected early-time occurred within less 
than two points from a zero crossing. If this observation 


proves to be generally true in later research, it may serve 


IEA 


as a way to check the starting point selected for one 


algorithm in terms of the other. 


C. CONCLUSIONS 

Both the Kumaresan-Tufts and the Cadzow-Solomon algorithms 
can effectively extract poles from the scattering response of 
a radar target. Because both algorithms obtain a least- 
squares solution to the system model, both perform acceptably 
in the presence of noise. Although eigenvalue compensation 
is not analytically justified in the Cadzow-Solomon algorithn, 
the results obtained through eigenvalue compensation in this 
method were generally superior to those similarly obtained in 
the Kumaresan-Tufts method. The results demonstrated the 
inherent advantages of an algorithm capable of processing a 
target's strongest response in the early time. The Kumaresan- 
Tufts method compared favorably with the Cadzow-Solomon only 


in responses with a long late—time. 


APPENDIX A. 


The following program implements’7 the 


THE KUMARESAN-TUFTS POLE EXTRACTION ALGORTITHM 


Kumaresan-Tufts 


algorithm as described in Chapter 2 of this thesis. The 


program is written in Fortran 77. 


The SVD and root-finding 


subroutines called by this program are found in the EISPACK 


itprary [18]. 


The SVD subroutine is a translation from ALGOL 


as given in [19]. The matrix multiplication and graphics 


subroutines, also called by this’ progran, 


Appendix C and D respectively. 


a 


INTEGER IERR,Kd,M,MN,MAGPOL, NSTRTPT , DELTAY 

INTEGER ITER ,NCAUS, NMENU,L/1/ 

INTEGER*2 KdPLT 

REAL*8 A(70,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1 (70) 

REAL*8 VS(70,70) ,UT(70,70) ,AINV(70,70) ,X(70) 

REAL*8 XP(70) ,B(70) ,SIGMA(70,70) ,SIG(70,70) 

REAL*8 COF (70) ,ROOTR(70) , ROOTI (70) 

REAL*8 D(1024) ,AVG,MACHEP/1.0E-16/ , Dy(140) 

COMPLEX*16 S(70) 

LOGICAL MATU/.TRUE./,MATV/.TRUE./ ,CAUSAL/.TRUE./,LONG/.TRUE. / 
LOGICAL DSET/. FALSE. /,NUFILE/ . TRUE. / 

CHARACTER TITLE*16 , HEADER*64 , YN*1, DC*1 , TITLER*16 , TITLEI*16 
CHARACTER TITL*16 


Fnter parameters for processing 


IF (DSET) CLOSE(10) 

NOVERLAY=0 

OPEN (10, FILE='PLOT') 

IF (DSET) GO TO 85 

WRITE (*,*) ‘Welcome to signal processing using the’ 
WRITE (*,*) ‘Kumaresan-Tufts method' 

WRITE (x, *) ' 8 

WRITE (*,*) ‘Do you want ' 

WRITE (x, *) ¢ J 

WRITE (*,*) '1. The long version for beginners' 
WRITE (*,*) '2. The short version for pros' 
WRITE (*,*) | ' 


pes 


are found in 


ibs) 


16 


10 


WRITE (*,*) ‘Please enter 1 or 2 ' 
READ (*,*) N 

IF (N .B). 1) THEN 

LONG=. TRUE. 

ELSEIF (N .EO. 2) THEN 
LONG=.FALSE. 

ELSE 

GO TO 15 

ENDIF 


WRITE (*,*) "Session will begin with entry of parameters needed fotr processing' 


WRITE (*,*) 

WRITE (*,*) ‘Do you want to enter parameters fron’ 
WRITE (*,*) q ] 

WRITE (*,*) '1. The keyboard’ 

WRITE (*,*) '2. A previously created file of parameters’ 
WRITE (*,*) ' ' 

WRITE (*,*) 'Please enter 1 or 2 ' 

READ (*,*) N 

IF (N .—EQ. 1) THEN 

GTO 1 

ELSEIF (N .BQ. 2) THEN 

WRITE (*,*) ‘Enter title of file containing parameters' 
READ (*,105) TITL 

OPEN (1, FILE=TITL) 

READ(1,105) TITLE 

READ (1,110) NPTS 

READ (1,110) NRT 

READ(1,110) Kd 

READ (1,110) M 

READ(1,110) DELTAY 

READ (1,110) NSTRTPT 

READ(1,110) NCAUS 

CLOSE (1) 

GO TO 85 

ELSE 

GO TO 16 

ENDIF 

WRITE (x, *) ' ¢ 


NUFILE=. TRUE. 

IF (.NOT. DSET) NSTRTPT=1 

WRITE (*,*) ‘Enter title of data file to be read’ 

READ (*,105) TITLE 

OPEN (1 , FILE=TITLE) 

READ (1,105) HEADER 

READ(1,110) NPTS 

IF (NPTS .GT. 1024) THEN 

WRITE (*,*) ‘Number of points in data file exceeds the dimension' 
WRITE (*,*) ‘of the array used in the program to store the file’ 
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20 


25 


STOP 
ENDIF 
CLOSE (1) 


IF (DSET) THEN 
IF (NSTRTPT+(Kd+M-1) *DELTAY .LE. NPTS) GO TO 85 
ENDIF 


IF (NUFILE) THEN 
WRITE (*,*) ‘Enter Kd, >= the estimated order of the systen ' 
READ (*,*) Kd 
IF (Kd .GT. 69) THEN 
WRITE (*,*) 'Kd must be less than 70, or dimension statements’ 
WRITE (*,*) ‘in this program must changed by the user' 
GO TO 3 
ELSEIF (Kd .LT. 2) THEN 
WRITE (*,*) ‘Kd must be at least 2' 
GO TO 3 
ENDIF 
IF (2*Kd .GT. NPTS) THEN 
WRITE (*,*) 'Kd must be less than or equal to ',NPTS/2 
GO TO 3 
ELSEIF (2*Kd .E). NPTS) THEN 
WRITE (*,*) ‘Kd equals' ,Kd 
WRITE (*,*) 'M must be',Kd 
M=Kd 
WRITE (*,*) ‘since there are a total of',NPTS 
WRITE (*,*) ‘points in ',TITLE 
GO TO 45 
ENDIF 
GO TO 4 
ELSEIF (DSET) THEN 
N=M 
IF (NSTRTPT+ (N+M-1)*DELTAY .LE. NPTS) THEN 
WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) 'Kd may range from ' NRT 
WRITE (*,*) ' to',N 
WRITE (*,*) ‘Enter Kd' 
READ (*,*) Kd 
IF (Kd .GE. NRT .AND. Kd .LE. N) GO TO 85 
GO TO 25 
ELSE 
NeN-1 
GO TO 20 
ENDIF 
ENDIF 


IF (NUFILE) THEN 


WRITE (*,*) ‘Enter M, the row dimension of the data matrix’ 
IF (.NOT. DSET .AND. LONG) THEN 


rr5 


30 


Sis 


45 


Sa 


WRITE (*,*) ' ' 

WRITE (*,*) ‘Note: Kd+M points in ',title 
WRITE (*,*) ' will be processed ' 
WRITE (x, *) * 

ENDIF 

WRITE (*,*) 'M may range from',Kd 

IF (NPTS-Kd .GT. 69) THEN 


WRITE (*,*) ' to 69' 
ELSE 

WRITE (*,*) ‘ to’ ,NPTS-Kd 
ENDIF 

READ (*,*) M 


IF (M .GT. 69) THEN 
WRITE (*,*) '"M must also be less than 70' 
GO TO 30 
ELSEIF (M .LT. Kd) THEN 
WRITE (*,*) 'M must be greater than or equal to Kd, Kd= ',Kd 
GO TO 30 
ELSEIF (Kd+4 .GT. NPTS) THEN 
WRITE (*,*) ‘"Kd+M must be less than or equal to',NPTS,',' 
WRITE (*,*) "the number of data points in',TITLE 
WRITE (*,*) ' ' 
G TO 30 
ENDIF 
ELSE 
N=Kd 
IF (NSTRTPT+ (Kd+N-1) *DELTAY .LE. NPTS) THEN 
NEN+1 
GO TO 35 
ELSE 
N=N-1 
ENDIF 
IF (N .BD. Kd) THEN 
WRITE (*,*) ‘M must equal',Kd 
M=Kd 
GO TO 85 
ENDIF 
IF (N .GT. 69) N-69 
WRITE (*,*) 'M may range from',Kd 
WRITE (*,*) ' to',N 
WRITE (*,*) ‘Enter M' 
READ (*,*) M 
IF (M .GE. Kd .AND. M .LE. N) GO TO 85 
G TO 40 
ENDIF 


IF (.NOT. NUFILE) GO TO 85 


N=1 
IF (NSTRTPT+N* (Kd+4-1) .LE. NPTS) THEN 


Lr6 


20 


N=N+1 

GO TO 50 

ELSE 

N=N-1 

ENDIF 

IF (N .BQ. 1) THEN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) "Spacing can only be 1' 

DELTAY=1 

IF (NUFILE) THEN 

GO TO 60 

ELSE 

GO TO 85 

ENDIF 

ENDIF 

IF (.NOT. DSET .AND. LONG) THEN 

WRITE (*,*) ‘Enter spacing between the ',Kd+¥ 
WRITE (*,*) ‘data points of ',TITLE 

WRITE (*,*) ‘to be processed ' 

WRITE (*,*) "4 

WRITE (*,*) ‘If, for example, one is chosen, then ',Kd+" 
WRITE (*,*) ‘consecutive points in ',TITLE 
WRITE (*,*) ‘will be processed ' 


WRITE (*, x) 

ENDIF 

WRITE (*,*) ‘Spacing may range fron 1' 
WRITE (*,*) ' to',N 


READ (*,*) DELTAY 

IF (DELTAY .GE. 1 .AND. DELTAY .LE. N) THEN 
IF (NUFILE) THEN 

GO TO 60 

ELSE 

GO TO 85 

ENDIF 

ELSE 

GO TO 55 

ENDIF 


WRITE (*,*) 'Do you wish to adjust eigenvalues? (y/n)' 
READ (*,120) YN 

IF (YN .&Q. 'N' .OR. YN .BD. 'n') THEN 

IF (NUFILE) GO TO 6 

GO TO 85 

ENDIF 

IF (YN .NE. 'Y' .AND. YN .NE. 'y') GO TO 60 

WRITE (*,*) 'Discard or compensate eigenvalues? (d/c)' 
READ (*,120) DC 

IF (DC .B. 'D' .OR. DC .&. 'd') G T 65 

IF (DC .NE. 'C' .AND. DC .NE. 'c') GOTO 2 

WRITE (*,*) ‘Enter estimate of the actual order of the system’ 


id 


65 


15 


WRITE (*,*) * ' 
IF (LONG) THEN 
WRITE (*,*) 'This estimate will be used to determine the ' 
WRITE (*,*) ‘number of eigenvalues compensated or discarded ' 
ENDIF 


WRITE (*,*) ‘the estimate may range from 2: 
WRITE (*,*) ' to’ ,Kd-1 
READ (*,*) NRT 

IF (NRT .GT. Kd .OR. NRT .LT. 2) THEN 

GO TO 65 

ELSEIF (.NOT. NUFILE) THEN 

GO TO 85 

ENDIF 

NSTRTPT=1 

IF (NSTRTPT+(Kd+-1)*DELTAY .LE. NPTS) THEN 
NSTRTPT=NSTRIPT+1 

G TO 70 

ELSE 

NSTRTPT=NSTRTIPT-1 

ENDIF 


IF (NSTRTPT .BQ. 1) THEN 

WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) ‘the starting point for processing the data’ 
WRITE (*,*) ‘must be the first point in the data file’ 
GO TO 85 

ENDIF 

WRITE (*,*) ‘Enter desired starting point in data file' 
IF (.NOT. DSET .AND. LONG) THEN 

WRITE (*,*) ‘1 indicates the first point in the data file ' 
ENDIF 

WRITE (* *) vas 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) ‘the starting point may range from die 
WRITE (*,*) ' to’ ,NSTRTPT 
READ (*,*) N 

IF (N .GE. 1 .AND. N .LE. NSTRTPT) THEN 

NSTRTPT=N 

ELSE 

WRITE (*,*) ‘Enter starting point again’ 

WRITE (*,*) * ' 

GO TO 75 

ENDIF 

IF (.NOT. NUFILE) GO TO 85 


WRITE (*,*) ‘Do you want the data matrix arrangement to be' 
WRITE (x, *) , 4 

WRITE (*,*) ‘1. Causal' 

WRITE (*,*) '2. Non-causal' 

WRITE (*,*) we 
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80 


12 


85 


WRITE (*,*) 'Please enter 1 or 2 ' 
READ (*,*) NCAUS 

IF (NCAUS .BD. 1) THEN 

CAUSAL=. TRUE. 

ELSEIF (NCAUS .BQ. 2) THEN 


G TO 85 

WRITE (*,*) ‘Enter title of file to contain parameters’ 
READ (*,105) TITL 
OPEN (1, FILE=TITL) 
WRITE (1,105) TITLE 
WRITE(1,110) NPTS 
WRITE (1,110) NRT 
WRITE(1,110) Kd 
WRITE(1,110) M 
WRITE(1,110) DELTAY 
WRITE (1,110) NSTRTPT 
WRITE (1,110) NCAUS 
CLOSE (1) 

IF (DSET) GO TO 85 


IF (DSET) THEN 

CLOSE (2) 

CLOSE (3) 

CALL SUBPLT (NOVERLAY) 
ENDIF 


DSET=.TRUE. 

NUFILE=.FALSE. 

WRITE (*,*) ' ' 

WRITE (*,*) '1. Data file to be processed or 
+ITLE 

WRITE(*,*) ' Number of data points in data file ' NPTS 
WRITE (*,*) '2. Estimated order of the system ' NRT 
WRITE (*,*) '3. Kd, the number of columns in the data matrix',Kd 
WRITE (*,*) '4. M, the number of rows in the data matrix',M 
WRITE(*,*) '5. Spacing between data points being processed ',DELTA 
+Y 

WRITE(*,*) '6. First point in the data file to be processed' ,NSTRT 
+PT 


WRITE (*,*) ' Last point in the data file to be processed’ ,NSTRT 
+PT+Kd+M-1 

IF (NCAUS .BQ. 1) THEN 

WRITE(*,*) '7. Data matrix arrangement for processing CA 
+USAL 

ELSE 
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90 


95 


100 


105 
110 
115 


WRITE(*,*) '7. Data matrix arrangement for processing NON-CA 
+USALsS' 

ENDIF 

WRITE (*, *) Ss] ' 

WRITE(*,*) '8. Begin processing using above settings’ 

WRITE (*,*) ‘9. Store parameters 1-7 in a file' 

WRITE (*,*) ‘10. Retrieve parameters 1-7 from a previously created 
+file' 

WRITE (*,*) ‘'11. Reset overlays' 

WRITE (*,*) '12. Re-plot overlays’ 

WRITE (*,*) ‘13. End this session of Kumaresan-Tufts signal process 
+ing' 

WRITE(*,*) ° ' 

WRITE (*,*) ‘Enter an integer from 1 to 12 to make changes as often 
+ as you desire’ 

READ (*,*) NMENU 

IF (NMENU .LT. 1 .OR. NMENU .GT. 13) THEN 

WRITE(*,*) ‘Enter an integer from 1 to 13' 

G TO 90 

ENDIF 


G TO (1,2,3,4,5,6,7,8,9,10,11,12,13) ,NMENU 


OPEN (1, FILE=TITLE) 
READ(1,105) HEADER 
READ(1,110) NPTS 
READ (1,115) XO 
READ(1,115) XQ 

DO 95 I=1,NPTS 
READ (1,115) D(I) 
CONTINUE 

CLOSE (1) 

KdPLT=Kd 
WRITE(*,*) ‘Enter title of file to contain real part of poles' 
READ (*,105) TITLER 
OPEN (2, file=TITLER) 


WRITE (*,*) "Enter title of file to contain imaginary part of poles' 
READ (*,105) TITLET 

OPEN (3, £1 le=-TITLET) 

WRITE(10,100) (KdPLT) 

WRITE (10,105) TITLER 

WRITE (10,105) TITLET 

FORMAT (12) 


MN=MAX (M, Kd) 
FORMAT (A) 


FORMAT (15) 
FORMAT (E12. 6) 


120 


120 


125 


130 


135 
140 


145 


150 


155 


160 


165 


FORMAT (A1) 
Form data matrix 


DO 125 I=1,Kd+M 
Dy (I)=D( (I-1) *DELTAY+NSTRTPT) 
CONTINUE 


DO 140 I=1,M 
DO 135 J=1,Kd 
A(I,J)=Dy (I+J) 
CONTINUE 
CONTINUE 


B(1)=Dy (1) 

DO 145 I=2,M 
B(I)=A(I-1,1) 
CONTINUE 


Begin singular value decomposition 
CALL SVD (MACHEP,M,Kd,MN,A,W,MATU,U,MATV, V, TERR, RV1) 


Errors in SVD? 

IF (IERR .GT. 0.0) THEN 

WRITE (*,*) ‘Error in singular value number ', IERR, STOP 
ENDIF 

IF (YN .—Q. 'N') G TO 190 


DO 150 I=1,Kd 
XP (T)=0.0 
CONTINUE 


Discard or compensate eigenvalues 
Order singular values 


XP (1) =W (1) 

DO 165 I=2,Kd 

DO 160 J=1,I 

IF (W(I) .GT. XP(J)) THEN 
DO 155 K=I+1,J,-1 
XP (K) =XP (K-1) 

XP (J) =W (T) 

GO TO 165 

ENDIF 

CONTINUE 

XP (I+1) =W (I) 
CONTINUE 


XP( ) now contains ordered singular values-XP(1) is the largest 


Mya 


170 


175 


180 
185 


190 


195 


205 
210 


Discard eigenvalues 
IF (DC .BQ. 'D') THEN 
DO 170 J=NRT+1,Kd 
W(J)=(0.0) 

ELSE 

Compensate eigenvalues 
AVG=0.0 

DO 175 J=NRT+1,Kd 
AVG=AVGHXP (J) **2 
CONTINUE 

IF (Kd .GT. NRT) AVG=AVG/DBLE (FLOAT (Kd-NRT) ) 


DO 185 J=1,Kd 

DO 180 K=1,Kd 

IF ( W(J) .BQ. XP(K) ) THEN 
IF ( K .GT. NRT ) THEN 
W(J)=0.0 
ELSE 
W(J)=DSORT (DABS ( W (J) *W(J) -AVG) ) 
ENDIF 
GO TO 185 
ENDIF 

CONTINUE 

CONTINUE 

ENDIF 


DO 200 I=1,M 

DO 195 J=1,M 

UT (I,J)=(U(J,I)) 
CONTINUE 
CONTINUE 


Form SIGMA+ (KdxM) 
DO 210 I=1,Kd 

DO 205 J=1,M 

SIGMA (I,J)=0.0 

IF (I .BQ. J .AND. W(J) .NE. 0.0) THEN 
SIGMA (I,J) =1.0D0/W (J) 
ELSE 

SIGMA (I,J)=0.0d0 
ENDIF 

CONTINUE 

CONTINUE 


Form SIGMA (MxKd) 

DO 220 I=1,M 

DO 215 J=1,Kd 

SIG(I,J)=0.0 

IF (I .BQ. J) SIG(I,J)=W(J) 


£22 


ZA5 
220 


229 


230 


235 


CONTINUE 
CONTINUE 


V=Kdxkd , SIGMA+=KdxM , VS=Kdax 
CALL MXMUL(V, SIGMA,Kd,Kd,M, VS) 


VS=KdxM , UT=MxM, AINV=KdxM 
CALL MXMUL (VS ,UT,Kd,M,M,AINV) 


Calculate matrix multiplication of AINV x B, where 
AINV=KdxM , B=Mx1 , XP=Kdx1 
CALL MXMUL (AINV,B,Kd,M,L, XP) 


Calculate autoregressive coefficients from prediction coefficients 
IF (XP(Kd) .EQ. 0.0) THEN 
WRITE (*,*) ‘ERROR, avoiding division by zero' 


STOP 

ELSE 

B(Kd)=1.0d0/XP (Kd) 
ENDIF 

DO 225 I=2,Kd 

B (I-1) =-B (Kd) *XP (Kd-I+1) 
CONTINUE 

DO 230 I=1,Kd 


X(I)=-B (Kd-I+1) 

IF (NCAUS .BQ. 1) X(I)=-XP(Kd-I+1) 
CONTINUE 

X(Kd+1)=1.0 


Compute the roots of the polynomial in z 
CALL POLRT(X,COF,KD,ROOTR, ROOTI, IER) 
IF (IER .NE. 0) WRITE (*,*) ‘ERROR with POLRT, IER=',IER,STOP 


DO 235 I=1,Kd 

WRITE (2,115) ROOTR(I) 

WRITE (3,115) ROOTI(I) 

S (I) =DCMPLX (ROOTR (I) , ROOTI (I) ) 
CONTINUE 


MAGPOL=0 

DO 240 I=1,Kd 

IF (CDABS(S(I)) .GE. 1.0d0) MAGPOL=MAGPOL+1 
CONTINUE 


WRITE (*,*) '# of poles with magnitude <= 1',Kd-MAGPOL 


WRITE (*,*) ‘HIT ANY KEY TO CONTINUE’ 
READ (*,105) HEADER 


ns 


245 


13 


Plot poles 


NOVERLAY=NOVERLAY+1 
CLOSE (2) 

CLOSE (3) 

CALL SUBPLT (NOVERLAY) 


J=0 
K=0 


DO 245 I=1,Kd 

IF (CDABS(S(I)) .LT. 1.0) THEN 
J=J+1 

K=K+1 

WRITE (*,*) S(I) ,CDABS(S(I)) 
ENDIF 

IF (J .—Q. 20) THEN 

WRITE (*,*) ‘Enter any key to continue' 
READ (*,105) HEADER 

J=0 

ENDIF 

CONTINUE 


WRITE(*,*) ‘Poles with magnitude less than one: ',K 
GO TO 85 


STOP 
END 
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APPENDIX B: THE CADZOW-SOLOMON POLE EXTRACTION ALGORITHM 


The following program implements the Cadzow-Solomon 
algorithm as described in Chapter 2 of this thesis. The 
program is written in Fortran 77. The SVD and root-finding 
subroutines called by this program are found in the EISPACK 
library {18]. The SVD subroutine is a translation from ALGOL 
as given in [{19]. The matrix multiplication and graphics 
subroutines, also called by this program, are found in 


Appendix C and D respectively. 


SLARGE 
INTEGER TERR, Kd, Kn,M,MN,MAGPOL, NSTRTPT , DELTAY 
INTEGER IER,NCAUS ,NMENU, INSTRTPT 
INTHEGER*2 KGPLT 
REAL*8 A(70,70) ,W(70) ,U(70,70) ,V(70, 70) ,RV1(70) 
REAL*8 VS(70,70) ,UT(70,70) ,AINV(70,70) ,X(70) 
REAL*8 XP (70) ,B(70) , SIGMA(70,70) ,SIG(70, 70) 
REAL*8 COF (70) ,ROOTR(70) ,ROOTI (70) 
REAL MAG 
REAL*8 D(1024) ,AVG, MACHEP/1.0E-16/ , Dy (140) ,Dx (1024) 
COMPLEX*16 S(70) 
LOGICAL MATU/ .TRUE./,MATV/.TRUE./,CAUSAL/. TRUE. /, LONG/. TRUE. / 
LOGICAL DSET/.FALSE. / ,NUFILE/ .TRUE. / 
CHARACTER TITLE*16, HEADER*64 , YN*1 ,DC*1, TITLER*16 , TITLEI*16 
CHARACTER TITL*16 , TITLD*16 


C Enter parameters for processing 


14 IF (DSET) CLOSE(10) 
NOVERLAY=0 
OPEN (10, FILE='PLOT' ) 
IF (DSET) GO TO 215 
WRITE (*,*) ‘Welcome to signal processing using the’ 
WRITE (*,*) 'Cadzow-Soloamon method' 
WRITE Gan*) 7 ' 
WRITE (*,*) 'Do you want ' 


lez 


aS 


35 


13 


WRITE (x *) 1 4 

WRITE (*,*) ‘1. The long version for beginners’ 
WRITE (*,*) '2. The short version for pros' 
WRITE (*, *) ve 

WRITE (*,*) ‘Please enter 1 or 2 ' 

READ (*,*) N 

IF (N .BQ. 1) THEN 

LONG=. TRUE. 

ELSEIF (N .BQ. 2) THEN 

LONG=.FALSE. 

ELSE 

GO TO 25 

ENDIF 


WRITE (*,*) ‘Session will begin with entry of parameters needed fo 


+r processing' 


WRITE. (*,*) 
WRITE (*,*) "Do you want to enter parameters fron' 
WRITE (x *) t ' 

WRITE (*,*) '1. The keyboard' 

WRITE (*,*) '2. A previously created file of parameters' 
WRITE (*,*) ' 

WRITE (*,*) ‘Please enter 1 or 2 ' 

READ (*,*) N 

IF (N .—O. 1) THEN 

G TO 8 

ELSEIF (N .EQ. 2) THEN 

WRITE (*,*) 'Enter title of file containing parameters’ 
READ (*,100) TITL 

OPEN (1 , FILE=TITL) 

READ (1,100) TITLE 

READ (1,110) NPTS 

READ (1,110) NRT 

READ(1,110) Kd 

READ(1,110) M 

READ(1,110) DELTAY 

READ (1,110) NSTRTPT 

READ(1,110) NCAUS 

READ(1,100) TITLD 

READ (1,110) NDPTS 

READ(1,110) Kn 

READ (1,110) INSTRTPT 

CLOSE (1) 

GO TO 215 

ELSE 

GO TO 35 

ENDIF 

WRITE (*,*) ' ' 


WRITE (*,*) ‘Enter title of file containing excitation waveform' 


126 


45 


10 
55 


READ (*,100) TITLD 

OPEN (8 , FILE=TITLD) 

READ (8,100) HEADER 

READ (8,110) N 

IF (N .GT. 1024) THEN 

WRITE (*,*) 'Number of points in data file exceeds the dimension’ 
WRITE (*,*) ‘of the array used in the program to store the file' 
STOP 


ENDIF 

CLOSE (8) 

IF ((N .GE. NDPTS) .AND. DSET) THEN 
NDPTS=N 

GO TO 215 

ENDIF 

NDPTS=N 


WRITE (*,*) ‘Enter estimated order of waveform’ 

IF (DSET) THEN 

MAXIMUM=NDPTS-M 

IF (MAXIMUM .GT. M-Kd-1) MAXIMUM=*-Kd-1 

IF (MAXIMUM .GT. NDPTS-INSTRTPT-Kn-M+1) THEN 
MAXIMUM=NDPTS-INSTRTPT-Kn-M+1 

ENDIF 

ELSE 

MAXTMUM=66 

ENDIF 

IF (MAXIMUM .ED. 1) THEN 

WRITE (*,*) 'The estimated order of the waveform can only be 1' 
IF (DSET) GO TO 215 

GO TO 10 

ELSE 

IF (DSET) THEN 

WRITE (*,*) ‘Given the other parameters chosen thus far,' 
ENDIF 

WRITE (*,*) ‘the order may range from ie 
WRITE (*,*) ' to’ , MAXIMUM 
READ (*,*) Kn 

IF (Kn .GE. 1 .AND. Kn .LE. MAXIMUM) THEN 

IF (DSET) GO TO 215 

GO TO 10 

ENDIF 

WRITE (*,*) ‘Enter estimated order again’ 

WRITE (* x) ? ' 

GO TO 45 

ENDIF 

IF (DSET) GO TO 215 


INSTRTPT=1 


IF (INSTRIPT+Kn+4-1 .GT. NDPTS) THEN 
INSTRTPT=INSTRIPT-1 


AT 


65 


ELSE 
INSTRTPT=INSTRTPT+1 
GO TO 55 

ENDIF 
MSTRT=INSTRTPT 


IF (INSTRIPT .Q. 1) THEN 

WRITE (*,*) ‘The first point can only be 1' 

GO TO 215 

ELSE 

WRITE (*,*) ‘Enter first point in waveform file to be processed’ 
WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) ‘the starting point may range from ih 
WRITE (*,*) ' to’ ,MSTRT 
READ (*,*) INSTRTPT 

IF (INSTRTPT .GE. 1 .AND. INSTRTPT .LE. MSTRT) THEN 

IF (DSET) GO TO 215 

GT 1 

ENDIF 

WRITE (*,*) ‘Enter starting point again' 

WRITE (*,*) " ' 

G TO 65 

ENDIF 

IF (DSET) GO TO 215 


IF (.NOT. DSET) NUFILE=.TRUE. 

IF (.NOT. DSET) NSTRTPT=1 

WRITE (*,*) ‘Enter title of data file to be read' 

READ (*,100) TITLE 

OPEN (12 , FILE=TITLE) 

READ (12,100) HEADER 

READ (12,110) NPTS 

IF (NPTS .GT. 1024) THEN 

WRITE (*,*) ‘Number of points in data file exceeds the dimension’ 
WRITE (*,*) ‘of the array used in the program to store the file’ 
STOP 

ENDIF 

CLOSE (12) 


IF (NUFILE) THEN 

GO TO 3 

ELSEIF (NSTRTPT+(Kd+M-1)*DELTAY .LE. NPTS) THEN 
GO TO 215 

ELSE 

GO TO 6 

ENDIF 


IF (NUFILE) THEN 
MAXIMUM=69-Kn-1 


28 


85 


25 


IF (MAXIMUM .GT. NPTS-69) MAXIMUM=NPTS-69 
IF (MIN .5Q. MAXIMUM) THEN 


WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) 'Kd must be ',MIN 


WRITE (*,*) ‘Enter Kd, >= the estimated order of the system ' 


WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) 'Kd may range from',MIN 

WRITE (*,*) ' to' , MAXIMUM 

READ (*,*) Kd 

IF (Kd .GE. MIN .AND. Kd .LE. MAXIMUM) @ TO 4 

GO TO 75 


ELSEIF (DSET) THEN 

MAXIMUM=M-Kn-1 

IF (MAXIMUM .GT. NPTS-M) MAXIMUM=NPTS-M 

MIN=2 

N=MAXIMUM 

IF (NSTRTPT+ (N+M-1) *DELTAY .LE. NPTS) THEN 

MAXIMUMEN 

IF (MIN .BQ. MAXIMUM) THEN 

Kd=MIN 

G TO 215 

ELSEIF (MAXIMUM .LT. MIN) THEN 

DELTAY=1 

IF (1+(2+4-1)*DELTAY .LE. NPTS) THEN 

Kd=2 

G TO 135 

ENDIF 

WRITE (*,*) ‘Error. Kd must be less than 2' 

Kd=2 

GO TO 215 

ENDIF 
WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) 'Kd may range from ' MIN 
WRITE (*,*) ' to’ ,MAXIMUM 
WRITE (*,*) ‘Enter Kd' 
READ (*,*) Kd 
IF (Kd .GE. MIN .AND. Kd .LE. MAXIMUM) GO TO 215 
GO TO 95 

ELSE 
NEN-1 
GO TO 85 

ENDIF 

ENDIF 


MS; 


C Determine M 

4 IF (NUFILE) THEN 
WRITE (*,*) 'Enter M, the row dimension of the data matrix’ 
IF (.NOT. DSET .AND. LONG) THEN 


WRITE (*,*) ° ' 

WRITE (*,*) 'Note: Kd+M points in ',title 
Wrae (70) will be processed ' 
WRITE (*,*) ' ' 

ENDIF 


105 WRITE (*,*) 'M may range from',Kd 
IF (NPTS-Kd .GT. 69) THEN 


WRITE (*,*) ' to 69' 
ELSE 

WRITE (*,*) ' to' ,NPTS-Kd 
ENDIF 

READ (*,*) M 


IF (M .GT. 69) THEN 
WRITE (*,*) 'M must also be less than 70' 
GO TO 105 
ELSEIF (M .LT. Kd) THEN 
WRITE (*,*) 'M must be greater than or equal to Kd, Kd= ',Kd 
GO TO 105 
ELSEIF (Kd+M .GT. NPTS) THEN 
WRITE (*,*) 'Kd+M must be less than or equal to',NPTS,',' 
WRITE (*,*) ‘the number of data points in',TITLE 
WRITE (x, *)  ] } 
GO TO 105 
ENDIF 
C Begin part for data already set 
ELSE 
N=Kd 
115 IF (NSTRTPT+ (Kd+N-1) *DELTAY .LE. NPTS) THEN 
NEN+1 
G TO 115 
ELSE 
NeN-1 
ENDIF 
IF (N .B). Kd) THEN 
WRITE (*,*) 'M must equal',Kd 
M=Kd 
GO TO 215 
ENDIF 
MAXIMUMEN 
IF (MAXIMUM .GT. 69) MAXIMUM69 
IF (Kd+Knt+1 .B). MAXIMUM) THEN 
M=Kd+Knt+1 
GO TO 215 
ELSEIF (Kd+Kn+1 .GT. MAXIMUM) THEN 
WRITE (*,*) ‘Kd must be reduced' 


ps0, 


125 


135 


145 


155 


G TO 3 

ELSE 

MINEKd+Knt+1 

ENDIF 
IF (MIN .LT. KntKd+]) MINeXKntKd+1 
WRITE (*,*) 'M may range from' ,MIN 
WRITE (*,*) ' to’ , MAXIMUM 
WRITE (*,*) ‘Enter M' 
READ (*,*) M 
IF (M .GE. MIN .AND. M .LE. MAXIMUM) GO TO 215 
G TO 125 
ENDIF 


Determine DELTAY 

IF (.NOT. NUFILE) GO TO 215 

h=1 

IF (NSTRTPT+N* (Kd+M-1) .LE. NPTS) THEN 
NeN+1 
GO TO 145 

ELSE 
N=N-1 

ENDIF 

IF (N .BD. 1) THEN 
WRITE (*,*) ‘Given the other parameters chosen thus far,' 
WRITE (*,*) ‘Spacing can only be 1' 
DELTAY=1 

IF (NUFILE) THEN 

GO TO 165 


IF (.NOT. DSET .AND. LONG) THEN 

WRITE (*,*) ‘Enter spacing between the ',Kd+¥ 

WRITE (*,*) ‘data points of ',TITLE 

WRITE (*,*) 'to be processed ' 

WRITE (*,*) ' ' 

WRITE (*,*) ‘If, for example, one is chosen, then ',Kd+4 
WRITE (*,*) ‘consecutive points in ',TITLE 

WRITE (*,*) ‘will be processed ' 

WRITE (*,*) ' ' 


WRITE (*,*) 'Enter spacing ' 

WRITE (*, *) ' t 

ENDIF 

WRITE (*,*) ‘Spacing may range fron be 
WRITE (*,*) ' to',N 

READ (*,*) DELTAY 

IF (DELTAY .GE. 1 .AND. DELTAY .LE. N) THEN 

IF (NUFILE) THEN 


de 


165 


175 


185 


GO TO 165 
ELSE 
G TO 215 
ENDIF 
ELSE 
GO TO 155 
ENDIF 


WRITE (*,*) ‘Do you wish to adjust eigenvalues? (y/n)' 

READ (*,150) YN 

IF (YN .—O. 'N' .OR. YN .BD. 'n') THEN 

IF (NUFILE) GO TO 6 

GO TO 215 

ENDIF 

IF (YN .NE. 'Y' .AND. YN .NE. ‘y') GO TO 165 

WRITE (*,*) ‘Discard or compensate eigenvalues? (d/c)' 

READ (*,150) DC 

IF (DC .FO. 'D' .OR. DC .B). 'd') THEN 

NRT=Kd 

GO TO 175 

ENDIF 

IF (DC .NE. 'C' .AND. DC .NE. 'c') GOTO 2 

WRITE (*,*) ‘Enter estimate of the actual order of the system’ 
WRITE (x) ? t 

IF (LONG) THEN 

WRITE (*,*) ‘This estimate will be used to determine the ' 
WRITE (*,*) ‘number of eigenvalues compensated or discarded ' 
ENDIF 


WRITE (*,*) 'the estimate may range fran Zi 
WRITE (*,*) ' to’ ,Kd+Knt1 
READ (*,*) NRT 

IF (NRT .GT. Kd+Kntl .OR. NRT .LT. 2) THEN 

GO TO 175 

ELSEIF (.NOT. NUFILE) THEN 

G TO 215 

ENDIF 

NSTRTPT=1 

IF (NSTRTPT+ (Kd+4-1) *DELTAY .LE. NPTS) THEN 
NSTRTPT=NSTRTIPT+1 

GO TO 185 

ELSE 

NSTRTPT=NSTRTIPT-1 

ENDIF 


IF (NSTRTPT .ED. 1) THEN 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) ‘the starting point for processing the data’ 
WRITE (*,*) ‘must be the first point in the data file’ 
G TO 215 

ENDIF 


iS 2 


iis 


205 


a2 


WRITE (*,*) ‘Enter desired starting point in data file' 
IF (.NOT. DSET .AND. LONG) THEN 

WRITE (*,*) '1 indicates the first point in the data file ' 
ENDIF 

WRITE (x ,*) i] i} 

WRITE (*,*) 'Given the other parameters chosen thus far,' 
WRITE (*,*) ‘the starting point may range from HE 
WRITE (*,*) ' to' ,NSTRTPT 
READ (*,*) N 

IF (N .GE. 1 .AND. N .LE. NSTRTPT) THEN 

NSTRTPT=N 

ELSE 

WRITE (*,*) ‘Enter starting point again' 

WRITE (* ,*) t ? 

GO TO 195 

ENDIF 

IF (.NOT. NUFILE) GO TO 215 


IF (DSET) THEN 

IF (NCAUS .ED. 1) THEN 
NCAUS=2 

GO TO 215 

ELSE 

NCAUS=1 

GO TO 215 

ENDIF 

ENDIF 

WRITE (*,*) 'Do you want the data matrix arrangement to be' 
WRITE (*,*) ' ' 

WRITE (*,*) '1. Causal’ 
WRITE (*,*) '2. Non-causal' 
WRITE (*%,*) ° ' 

WRITE (*,*) ‘Please enter 1 or 2 ' 
READ (*,*) NCAUS 

IF (NCAUS .BD. 1) THEN 
CAUSAL=. TRUE. 

ELSEIF (NCAUS .ED. 2) THEN 
CAUSAL=.FALSE. 

ELSE 

GO TO 205 

ENDIF 

GO TO 215 


WRITE (*,*) ‘Enter title of file to contain parameters' 
READ (*,100) TITL 

OPEN (1, FILE=TITL) 

WRITE (1,100) TITLE 

WRITE (1,110) NPTS 

WRITE (1,110) NRT 


33 


15 


Z¥5 


WRITE (1,110) Kd 
WRITE(1,110) M 
WRITE(1,110) DELTAY 
WRITE(1,110) NSTRTPT 
WRITE (1,110) NCAUS 
WRITE(1,100) TITLD 
WRITE (1,110) NDPTS 
WRITE(1,110) Kn 

WRITE (1,110) INSTRTPT 


CLOSE (1) 

IF (DSET) GO TO 215 

IF (DSET) THEN 

CLOSE (2) 

CLOSE (3) 

CALL SUBPLT (NOVERLAY) 

ENDIF 

DSET=.TRUE. 

NUFILE=.FALSE. 

WRITE(*,*) * * 

WRITE(*,*) ‘1. Data file to be processed ran 
+ITLE 

WRITE(*,*) ° Number of data points in data file ' NPTS 
WRITE(*,*) ‘2. Estimated order of the system ’ NRT 
WRITE (*,*) '3. Kd, the number of columns in the data matrix',Kd 
WRITE (*,*) '4. M, the number of rows’ in the data matrix',M 
WRITE (*,*) '5. Spacing between data points being processed ',DELTA 
+Y 

WRITE (*,*) '6. First point in the data file to be processed' ,NSTRT 
+PT 

WRITE(*,*) ' Last point in the data file to be processed’ ,NSTRT 
+PT+Kd+i-1 

IF (NCAUS .EO. 1) THEN 

WRITE (*,*) '7. Data matrix arrangement for processing CA 
+USAL 

ELSE 

WRITE(*,*) '7. Data matrix arrangement for processing NON-CA 
+USALS' 

ENDIF 

WRITE (* , *) A] A] 

WRITE(*,*) '8. File containing excitation waveform uA 
+ITLD 

WRITE(*,*) ' Number of data points in above file ' NDPTS 
WRITE (*,*) '9. Estimated order of the waveform ' Kn 
WRITE(*,*) '10. First point in the file to be ' 

WRITE (*,*) ' input into the data matrix " INSTR 
+TPT 


134 


229 


11 


Z05 


245 


WRITES 3) = 


WRITE(*,*) '11. Begin processing using above settings' 
WRITE(*,*) '12. Store parameters 1-10 in a file' 
WRITE (*,*) '13. Retrieve parameters 1-10 from a previously created 


+ file’ 


WRITE(*,*) ‘14. Reset overlays’ 
WRITE (*,*) ‘15. Re-plot overlays' 
WRITE(*,*) '16. End this session of Cadzow-Solomon signal processi 


+ng 


WRITE(*,*) ° ' 
WRITE(*,*) ‘Enter an integer from 1 to 16 to make changes as often 


+ as you desire' 


READ (*,*) NMENU 

IF (NMENU .LT. 1 .OR. NMENU .GT. 16) THEN 
WRITE (*,*) ‘Enter an integer from 1 to 16' 
GO TO 225 

ENDIF 


GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16) ,NMENU 


OPEN (12 , FILE=TITLE) 
READ(12,100) HEADER 
READ(12,110) NPTS 
READ (12,120) XQ 
READ (12,120) XQ 

DO 235 I=1,NPTS 
READ(12,120) D(I) 
CONTINUE 

CLOSE (12) 


OPEN (8 , FILE=TITLD) 
READ (8,100) HEADER 
READ (8,110) NDPTS 
READ (8,120) XQ 
READ (8,120) XQ 

DO 245 I=1,NDPTS 
READ(8,120) Dx(T) 
CONTINUE 

CLOSE (8) 


KdPLT=Kd 

WRITE (*,*) ‘enter title of file to contain real part of poles' 
READ (*, 100) TITLER 

OPEN (2, FILE=TITLER) 


WRITE (*,*) ‘enter title of file to contain imaginary part of poles’ 
READ(*,100) TITLEI 
OPEN (3 , FILE=TITLET) 


bo 


130 


100 
110 


150 


205 


265 


215 


285 


295 


305 


WRITE(10,130) (KdPLT) 
WRITE (10,100) TITLER 
WRITE(10,100) TITLEI 
FORMAT (12) 


MN=MAX (M, Kd+Knr+1) 


FORMAT (A) 
FORMAT (15) 
FORMAT (E12. 6) 
FORMAT (A) 


DO 255 I=1,Kd+M 
Dy (I) =D( (I-1) *DELTAY+NSTRTPT) 
CONTINUE 


DO 285 I=1,M 

DO 275 J=1,Kdtkn+1 

A(I,J)=Dy (I+J) 

IF (J .GE. Kdt1) A(I,J)=Dx (I+J+INSTRIPT-2-Kd) 
CONTINUE 

CONTINUE 


B(1) =Dy (1) 
DO 295 I=2,M 
B(I)=A(I-1,1) 
CONTINUE 


N=Kd+Kn+1 

Begin singular value decomposition 

CALL SVD (MACHEP,M,N,MN,A,W,MATU,U,MATV,V, IERR,RV1) 
Errors in SVD? 

IF (IERR .GT. 0.0) THEN 

WRITE (*,*) ‘Error in singular value number ', JERR, STOP 
ENDIF 

IF (YN .—D. 'N') GO TO 385 

DO 305 I=1,Kd+Kntl 

XP (I)=0.0 

CONTINUE 


Discard or compensate eigenvalues 
Order singular values 


XP (1) =W (1) 


136 


g15 


325 


335 


345 


855 


365 
B75 


385 


395 
405 


DO 335 I=2,Kd+Knt+l 
DO 325 J=1,I 

if (W(I) .GT. XP(J)) THEN 
DO 315 K=I+1,J,-1 
XP (K) =XP (K-1) 

XP (3) =W(1) 

GO TO 335 

ENDIF 

CONTINUE 

XP (I+1) =W(I) 
CONTINUE 


XP( ) now contains ordered singular values: XP(1) is the largest 


Discard eigenvalues 
IF (DC .BQ. 'D') THEN 
DO 345 J=NRT+1,Kd+Knt1 
W(J)=(0.0) 

ELSE 

Compensate eigenvalues 
AVG=0.0 

DO 355 J=NRT+1,Kd+Knt1 
AVG=AVG+XP (J) **2 
CONTINUE 

IF (Kd+Kn+] .GT. NRT) AVG=AVG/DBLE (FLOAT (Kd+Kn+1-NRT) ) 


DO 375 J=1,Kd+Kn+1 

DO 365 K=1,Kd+Knt1 

IF ( W(J) .BQ. XP(K) ) THEN 
IF ( K .GT. NRT ) THEN 
W(J)=0.0 
ELSE 
W(J)=DSORT (DABS ( W(J) *W(J) -AVG) ) 
ENDIF 
GO TO 375 
ENDIF 

CONTINUE 

CONTINUE 

ENDIF 


DO 405 I=1,M 

DO 395 J=1,M 
UT(I,J)=(U(J,I) ) 
CONTINUE 
CONTINUE 


Form SIGMA+ (Kd+Knt+1 x M) 


DO 425 I=1,Kd+Kntl 
DO 415 J=1,M 


UC, 


415 


435 
445 


455 


465 


SIGMA (I, J)=0.0 

IF (I .—Q. J .AND. W(J) .NE. 0.0) THEN 
SIGMA (I, J)=1.0d0/W (J) 

ELSE 

SIGMA (I,J)=0.0D0 

ENDIF 

CONTINUE 

CONTINUE 


Form SIGMA (M x KdtKn+1) 
DO 445 I=1,M 
DO 435 J=1,KdtKnt1 
SIG(I,J)=0.0 

IF (I .B). J) SIG(I,J)=W(J) 
CONTINUE 

CONTINUE 


V=Kd+Kn+1xkd+kn+1 , SIGMA+=Kd+kKn+ Lx, VS=Kd+Kn+LM 
CALL MXMUL(V, SIGMA, Kd+Knt+1 ,Kd+Kn+1,M, VS) 


VS=Kd+Knt+1eM , UT=MxM , AINV=Kd+Kn+1xM 
CALL MXMUL(VS,UT,Kd+Knt1 ,M,M, AINV) 


Calculate matrix multiplication of AINV x B, where 
AINV=Kd+Knt+LxM , B=Mxl , XP=Kd+Knt+1x1 
CALL MXMUL (AINV,B,Kd+Knt+1 ,M,L, XP) 


Compute autoregressive coefficients from prediction coefficients 


IF (XP(Kd) .BQ. 0.0) THEN 
WRITE (*,*) ‘ERROR, avoiding division by zero' 


STOP 

ELSE 

B (Kd) =1.0d0/XP (Kd) 
ENDIF 

DO 455 I=2,Kd 

B (I-1) =-B (Kd) *XP (Kd-1+1) 
CONTINUE 

DO 465 1=1,Kd 


X(I)=-B (Kd-I+1) 

IF (NCAUS .BQ. 1) X(I)=-XP (Kd-I+1) 
CONTINUE 

X (Kd+1)=1.0 

Compute the roots of the polynamial in z 
CALL POLRT (X,COF,KD,ROOTR , ROOTT , IER) 


IF (IER .NE. 0) WRITE (*,*) ‘ERROR with POLRT, IER=', IER, STOP 


ss 


475 


485 


495 


DO 475 I=1,Kd 

WRITE (2,120) ROOTR (I) 

WRITE (3,120) ROOTI(T) 

S (I) =DCMPLX (ROOTR (I) , ROOTI (I) ) 
CONTINUE 


MAGPOL=0 

DO 485 I=1,Kd 

IF (CDABS(S(I)) .GE. 1.0D0) MAGPOL=MAGPOL+1 
CONTINUE 


WRITE (*,*) '# of poles with magnitude <= 1',Kd-MAGPOL 
WRITE (*,*) ‘HIT ANY KEY TO CONTINUE’ 
READ (*,100) HEADER 


Plot poles 
NOVERLAY=NOVERLAY+1 
CLOSE (2) 

CLOSE (3) 

CALL SUBPLT (NOVERLAY) 


J=0 
K=0 


DO 495 I=1,Kd 

IF (CDABS(S(I)) .LT. 1.0) THEN 
WRITE (*,*) S(I) ,CDABS(S(I)) 
J=J+1 

K=K+1 

ENDIF 

IF (J .—Q. 20) THEN 

WRITE (*,*) "HIT ANY KEY TO CONTINUE' 
READ (*,100) HEADER 

Jz0 

ENDIF 

CONTINUE 


WRITE (*,*) ‘Poles with magnitude less than one ',K 
WRITE (*,*) ‘HIT ANY KEY TO CONTINUE’ 
READ (*,100) HEADER 


G TO 215 


END 


Tsay 
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APPENDIX C. MATRIX MULTIPLICATION 


SUBROUTINE MXMUL(A,B,RA,CA,CB, AB) 
INTEGER RA,CA,CB 
REAL*8 A(70,70) ,B(70,70) ,AB(70, 70) 


Calculates matrix multiplication of A x B=AB, where 
A=RAxCA, B=CAXCB , AB=RAXCB 


DO 30 I=1,RA 

DO 20 J=1,CB 

AB(I,J)=0.0 

DO 10 K=1,CA 
AB(I,J)=AB(I,J)+A(I,K) *B(K,J) 
CONTINUE 

CONTINUE 

CONTINUE 

RETURN 


END 
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APPENDIX D. GRAPHICS ROUTINE 


SUBROUTINE SUBPLT (NOVERLAY) 


MS-FORTRAN Program using "Grafmatic" Library Subroutines. 
Plots a Solid Line and Optional Overlay Plot for Comparison. 
Written by M.A. Morgan with Latest Update August 1989. 


Default Printer is "IBM Graphics" (e.g. Epson, Okidata, IBM) 
With Plot Rotated 90 degrees From the Vertical. "GrafPlus.Com" 
May be Run to Rotate Plot Upright on Paper and to Use a Variety 
of Impact Printers. "GrafLaser.Com'' May be Run to Use a Laser 
Printer. See GrafPlus/Laser Manual From Jewell Technology. 


CHARACTER*1 YN, YN1,DUM, YN2, SYMBOL, BELL, FEED, FFYN 
CHARACTER*4 LINE 

CHARACTER*7 SYMB 

CHARACTER*16 LTIT,CTIT, FNAME, TITLER , TITLE 
CHARACTER*64 TITLE, HCOPY 

REAL CRTR (70) ,CRTI (70) ,NRTR (70) ,NRTI (70) 
INTEGER*2 N,JROW, JCOL, ISYM1, ISYM2, ITYPE1 , ITYPE2 ,NSCRN 
INTEGER*2 CYAN, GREEN , WHITE, YELLOW, RED, BLACK, BLUE ,NTWO 
INTEGER*2 JROW1 , JROW2, JCOL1 , JCOL2, CROSS , KAPLT , I 
INTEGER*2 PURPLE, RUST 

EXTERNAL XFUN, YFUNP, YFUNN 

LINE='-— ) 

WHITE=7 

GREEN=10 

CYAN=11 

YELLOW=14 

RED=12 

BLACK=0 

BLUE=1 
NTWO=2 

PURPLE=5 

RUST=6 

BELL=CHAR (7) 

FEED=CHAR (12) 


Clear Screen and Put Up Introduction - on Blue Backgound for BGA 
Only; Another Background Color is Possible by Changing "BLUE" 
in the Calls to QPREG and QOVSCN. 

CALL OSMODE (NTWO) 

CALL QPREG(0,BLUE) 
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CALL QOVSCN (BLUE) 
WRITE(*,*) BELL 
NS=1 

NSCRN =16 
ITYPE2=0 


Calling GRAFMATIC Routines and Plotting Fi Solid Line Graph 
ITYPE1=0 
ISYM1=-1 
NDOTS1=0 
JROW1=1 
JROW2=350 
JCOLI= 75 
JCOL2= 565 
XMIN=-1 .2 
XMAX=1 .2 
YMIN=-1 .2 
YMAX=1 . 20 
YOVERX=1.115 
XORG=0 .0 
YORG=0 .0 
XST=-1.1 
XFIN=1.1 
YsTt—1.1 
YFIN=1.1 


CALL QSMODE (NSCRN) 
CALL QPLOT (JCOL1 , JCOL2, JROW1 , JROW2 , XMIN, XMAX, YMIN, YMAX, XORG, YORG, 


+1, YOVERX, 1.5) 


CALL QSETUP (NDOTS1, CYAN, ISYM1, RED) 


IF (XFIN-XST .LE. 9.0) XMAJOR=0.6 

IF (XFIN-XST .LE. 6.0) XMAJOR=0.4 

IF (XFIN-XST .LE. 3.3) XMAJOR=0.2 

IF (XFIN-XST .GE. 9.0) XMAJOR=(XFIN-XST) /10.0 
MINOR=0 

LABEL=1 

NDEC=2 

CALL QXAXIS (XST, XFIN, XMAJOR ,, MINOR , LABEL , NDEC) 
YMAJOR=XMAJOR 


CALL QYAXIS (YST, YFIN, YMAJOR , MINOR, LABEL , NDEC) 
Plot unit circle 

A=-1.0 

B=1.0 


CALL QCURV (XFUN, YFUNP A,B) 
CALL QCURV (XFUN, YFUNN, A,B) 


IF (NOVERLAY-1 .LT. 1) THEN 
IF (NOVERLAY-1 .BQ. 0) THEN 
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WRITE (*,3) NOVERLAY-1 
ELSE 

NZERO=0 

WRITE (*,3) NZERO 
ENDIF 

ELSEIF (NOVERLAY-1 .GT. 1) THEN 
WRITE (*,3) NOVERLAY-1 
ELSE 

WRITE (*,4) NOVERLAY-1 
ENDIF 

FORMAT (I3,' OVERLAYS’) 
FORMAT (I3,° OVERLAY ') 
REWIND (10) 


DO 20 I=1,NOVERLAY 

READ (10,110) KdPLT 

READ (10,100) TITLER 

READ (10,100) TITLEI 

OPEN (2 , FILE=TITLER) 

OPEN (3 , FILE=TITLEI) 

NKd=KdaPLT 

DO 27 J=1,KdPLT 

READ (2,120) NRTR(J) 

READ (3,120) NRTI (J) 

IF (DSORT (NRTR (J) **2+NRTI(J)**2) .GT. 1.1) THEN 
NKd=NKd-1 

NRTR (J) =0.0 

NRTT (J) =0.0 

ENDIF 

CONTINUE 

PURPLE=5 

RUST=6 

WHITE=7 

GREEN=10 

CYAN=11 

YELLOW=14 

RED=12 

BIUE=1 

IF (I .B). 1) THEN 

CALL QSETUP (NDOTS1 , CYAN, ISYM1 , RED) 
ELSEIF (I .BQ. 2) THEN 

CALL QSETUP (NDOTS1, CYAN, ISYM1 , GREEN) 
ELSEIF (I .BQ. 3) THEN 

CALL QSETUP (NDOTS1 , CYAN, ISYM1 , YELLOW) 
ELSEIF (I .BD. 4) THEN 

CALL QSETUP (NDOTS1,CYAN, ISYM1, BLUE) 
ELSEIF (I .BQ. 5) THEN 

CALL QSETUP (NDOTS1, CYAN, ISYM1 , WHITE) 
ELSEIF (I .BQ. 6) THEN 

CALL QSETUP (NDOTS1 , CYAN, ISYM1 , PURPLE) 
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ELSEIF (I .BQ. 7) THEN 


CALL QSETUP (NDOTS1 , CYAN, ISYM1 , RUST) 


ELSE 


CALL QSETUP (NDOTS1, CYAN, ISYM1,RED) 


ENDIF 
CALL QTABL (ITYPE1 , KAPLT,NRTR,NRTI) 
CONTINUE 


READ (*,100) DUM 

GO TO 40 

HCOPY='HARDCOPY——> ENTER P OR p' 
CALL QPTXT (30, COPY, RED, 25,1) 
CALL QCMOV (55,1) 

HCOPY=' 

CALL OPTXT (40, HCOPY, BLACK, 25,1) 


IF (DUM .NE. 'P' .AND. DUM .NE. 'p') GO TO 40 


CALL QPSCRN 

OPEN (1, FILE='PRN') 
WRITE (1,160) FEED 
FORMAT (A) 

FORMAT (12) 

FORMAT (E12.6) 
FORMAT(' ',A,\) 
CONTINUE 

CALL QSMODE (NTWO) 
CALL QPREG (0, BLUE) 
CALL QOVSCN (BLUE) 


WRITE (*,*) NKd, ‘points were plotted' 


RETURN 
END 


REAL FUNCTION XFUN(T) 
XFUN=T 

RETURN 

END 


REAL FUNCTION YFUNP (T) 
YFUNP=SORT (1 .0-T*T) 
RETURN 


END 


REAL FUNCTION YFUNN(T) 
YFUNN=-SORT (1 .0-T*T) 
RETURN 

END 
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